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1. The ideal Fermi gas 

We consider in this seetion a non- interacting Fermi gas in a single spin component, 
at thermal equilibrium in the grand canonical ensemble. We review basic properties of 
such a gas. We concentrate on the degenerate regime pX^ ^ 1, with p the density and A 
the thermal de Broglie wavelength defined as 



(1) = 



rnksT 



I'l. Basic facts. ~ In first quantization, the Hamiltonian of the non-interacting system 
is the sum of one-body terms, 



(2) H = J2hoi^) 



where N is the operator giving the number of particles, and ho is the one-body Hamil- 
tonian given here by 

(3) = 1^ + ^^""^ 

where p is the momentum operator, m is the mass of a particle and U{f) is a position 
dependent external potential. We shall consider two classes of external potential. Either 
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a zero external potential, supplemented by periodic boundary conditions in a box of 
size L, defined by f S [0,L['^, or a harmonic potential mimicking the traps used in 
experiments: 



(4) U{r) = J2 



r 

2 

Q = l 

where uja is the oscillation frequency of a particle in the trap along one of the d dimensions 
of space. 

The system is made of indistinguishable fermions, all in the same spin state. To imple- 
ment the antisymmetry of the many-body wavefunction, we move to second quantization. 
Introducing the orthonormal basis of the cigcnmodcs (px of the one-body Hamiltonian, 
with cigenenergy e^, we expand the field operator as 



(5) ,/'(f) = ^0A(r) 



where a\ annihilates one particle in the mode 0a and obeys anticommutation relations 
with the creation and annihilation operators: 

(6) {ax,ay) = Q 

(7) {ax,a\,} = 5x\'- 

The Hamiltonian then reads in second quantized form: 

(8) H = Y,^xAax. 

A 

and the operator giving the number of particles: 

(9) N = Y,a\ax. 

A 

The system is assumed to be at thermal equilibrium in the grand-canonical ensemble, 
with a many-body density operator given by 

(10) ^ ^ S-i 

where the factor S ensures that a has unit trace, (3 = l/ksT and /i is the chemical 
potential. The choice of the grand-canonical ensemble allows to decouple one mode 
from the other and makes the density operator Gaussian in the field variables, so that 
Wick's theorem may be used to calculate expectation values, which are thus (possibly 
complicated) functions of the only non-zero quadratic moments 



(11) {a\ax) = nx ^ n{ex) = 



exp[/3(eA - fJ-)] + l' 
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This is the famous Fermi-Dirac law for the occupation numbers n\. An elementary 
derivation of Wick's theorem can be found in an appendix of [1]. 

We note that, because of the fermionic nature of the system, the grand-canonical 
ensemble cannot be subject to unphysically large fluctuations of the number of particles, 
contrarily to the ideal Bose gas case. Using Wick's theorem we find that the variance of 
the particle number is 



(12) Var7V = ^nA(l 



which is indeed always below the Poisson value (N). In the zero temperature limit, for fi 
not coinciding with an eigenmode energy eA, one finds a vanishing variance of the particle 
number, so that the grand-canonical ensemble becomes equivalent to the canonical one. 

In this lecture, we shall be concerned with the degenerate limit, where the temperature 
is much smaller than the chemical potential: 

(13) fcsT<A«- 

The Fermi-Dirac distribution has then the shape presented in Fig. [1] We shall repeatedly 
use the fact that the occupation number of holes I —■ nx for e\ < ^ and the occupation 
number of particles n\ for e\ > fi are narrow functions of e\ — ii, of energy width ~ ksT. 

More precisely, a standard low-T expansion technique proceeds as follows [2]. One 
writes the Fermi-Dirac formula as 

/1.x / \ / N sign (e — u.) 

(14) n{e) = nT=o(e) + 



exp(/3|e-Ai|) + r 



singling out as a second term a narrow function of e — ^. In the case of a continuous 
density of states p(e), an integral over e appears in thermodynamic quantities. One then 
extends the integration over e to ] — oo, -|-oo[ for the second term of the right-hand side, 
neglecting a contribution 0[exp{— fj./ ksT)], after having expanded in powers of e — /i the 
functions of e multiplicating this second term, such as p{e). 

We illustrate this technique for the calculation of the density p in a spatially homo- 
geneous system in the thermodynamic limit: 



(15) Phom(Ai,T)= / TT^rifc = / dep{e)n{e) 



(27r)° 



+ 00 







where we introduced the free space density of states 



(16) p{e) = j j^S{e - h'e/2m) = -{2^)-''Va{V2^elh) 

with Vdiu) is the volume enclosed by the sphere of radius u in dimension d. We go to the 
limit T ^ for a fixed chemical potential, and we take /i > in order to have a non-zero 
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density at zero temperature. From the rewriting Eq. (jl4[) of the Fermi-Dirac formula we 
obtain 



(17) pho.(M, T) = p,^,, 0) + r diSe) + P^^' - 

Jo 



exp(/3(5e) + 1 



Oi 



where we cut the integration over e to 2/^, introducing an exponcntiaUy smaU error. The 
zero temperature contribution is readily evaluated by the integral in k space: 



(18) 



Phom(/i,0) = 



h^k'^/2m<fi 



d'^k _ Vrf(l) /2m/i\ 
(2^ ~ (2^ \W) 



d/2 



The thermal contribution is obtained as a low-T expansion by expanding p{fi ± 6e) in 
powers of Se: only odd powers of Se contribute; in the resulting integrals over Se, we 
extend the upper bound from p to +oo, paying the price of an exponentially small error. 
Finally a series expansion with even powers of T is obtained: 



(19) 



PhomifJ',T) = Phom(/i,0) 



O ((^bT/m)^) 



The bidimensional case d ~ 2 deserves a particular discussion. One finds that all the 
powers of T in the expansion have vanishing coefficients, since the density of states is 
constant, p{e) = ■m/{2'Kh'^). The integral over e in Eq. pS]) can be performed exactly, 



(20) 



ZTTh 



For a fixed and negative chemical potential, one sees that the density tends exponen- 
tially rapidly to zero when T ^ 0. For p > a more convenient rewriting in the low 
temperature limit is 



(21) 



Phom(M,r) 



mp 
2^ 



l + M^ln(l + e-^^) 
p 



This explicitly shows that the deviation of the spatial density from the zero temperature 
value Phom(Mi 0) is exponentially small for T ^ if d = 2. 

1'2. Coherence and correlation functions of the homogeneous gas. - We consider here 
a gas in a box with periodic boundary conditions. The eigenmodes are thus plane waves 
4'j:{r) = exp{ik ■ r)/L'^/^, with a wavevector equal to 2n/L times a vector with integer 
components. We shall immediately go to the thermodynamic limit, setting L to infinity 
for a fixed p and T. 

The first-order coherence function of the field is defined as the following thermal 
average: 



(22) 



giif) = {i^HrmO)). 
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In the thermodynamic hmit, it is given by the Fourier transform of the momentum 
distribution of the gas, which is an isotropic function: 

(23) g,ir ) = c^ir)^ J ^n^e^^^'^ 

It is easy to see that no fong-range order exists for this coherence function, even at 
zero temperature. At T = 0, rik = 1 for k < ko and rifc = for fc > fcg, where the 
wavevector ko is defined by 

and is simply the Fermi wavevector, in the present case of zero temperature. In the 
dimensions for d = 1 to d = 3 this leads to 

(25) '\)\D\r) = 

Trr 

(26) 02^(^) = ilLj^(fcor) 

zvrr 

(27) 03d(O = -^dAiD{r). 
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where Ji{x) is a Bessel function, behaving as — cos(a; + 7r/4)(2/7ra;)^/^ for x — > +00. One 
then gets an algebraic decay of (j> at large distances, in all dimensions, because of the 
discontinuity of the momentum distribution. Note that Eq. (|27p holds at all temperatures. 

At low but non-zero temperature, the Fermi surface is no longer infinitely sharp but 
has a smooth variation with an energy width ^ fcsT. This corresponds to a momentum 
width (5fco such that 

/oQ^ h^kpSko _ 

(28) e/co+Sfco - efco — = ks-L ■ 

m 

We then expect a decreasing envelope for the oscillating function <f>{r) with a decay 
faster than algebraic with a length scale 1/Sko. This can be checked by the following 
approximate calculation in ID: we restrict to /c > using parity, we single out from 
the zero-temperature contribution 6{ko — k), where 9 is the Heaviside function, and we 
linearize the kinetic energy dispersion relation around feg, writing CkQ+sk — iJ^ + hkoSk/m. 
Extending the integration over Sk to —00 we obtain 

29 (/.m r - (/.f^o r ~ sin fcoO / du \ °' 

TT Jo exp(M) -I- 1 

where u originates from the change of variable u — |5fc|/i5fco. Integrating by parts over u, 
taking the integral of the sine function, gives a fully integrated term exactly compensating 
the T = contribution, so that one is left with 

(30) <^,nir) ^ p^^ cos(ua-or) ^ 

^ ' ^ 27rr Jo cosh2(u/2) 

After extension to an integral over the whole real line (thanks to the parity of the inte- 
grand), one may close the integration contour at infinity and sum the residues over all 
the poles with Imu > to get: 

, , ^ Skp smjkor) 
smh(7roKor) 

We have recovered a result derived in |3] . This approximate formula can also be used to 
calculate the 3D case, by virtue of Eq. (|27|) . 

The correlation function of the gas, also called the pair distribution function, is defined 
as the following thermal average: 

(32) g2{r) - {tPHr)4'HomOWr))- 
By virtue of Wick's theorem, 52 is related to the function 4> as 

(33) ff2(r) = <^2(O)-0^(r) 
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where </>(0) is simply the mean particle density p. A key feature is that 172 vanishes in 
the limit r (quadratically in r): this is a direct consequence of the Pauli exclusion 
principle, and expresses the fact that one cannot find (with a finite probability density) 
two fermions with the same spin state at the same location. 

This spatial antibunching is intuitively expected to lead to reduced density fluctu- 
ations for an ideal degenerate Fermi gas, as compared to an ideal Bose gas. This is 
quantified in the next subsection. 

1'3. Fluctuations of the number of fermions in a given spatial zone. - We define the 
operator giving the number of particles within the sphere of radius R: 

(34) Nr= I d'^r{p\r){p{r), 

Jr<R 

the gas being spatially homogeneous and taken in the thermodynamic limit. We now 
show that the number N^. has subpoissonian fluctuations for a degenerate ideal Fermi 
gas. 

The mean value of the number of particles within the sphere is simply 

(35) {Nr) - pVd{K) 

where p = (f>{0) is the mean density and Vd{R) the volume enclosed by the sphere of 
radius R in dimension d. From Wick's theorem we find a variance 

(36) Var Nr={Nr)~[ d'^rj dV02(|f_f'|). 

Jr<R Jr'<R 

This immediately reveals the subpoissonian nature of the fluctuations. 

Wc quantify this subpoissonian nature in the limit of a large radius R. Replacing (jj 
by its expression as a Fourier transform of Uk, see Ea. (|23p . we obtain 

(37) YarNR={NR)^J ^ J n,n„ F {k ^ k') 

where F is the modulus squared of the Fourier transform of the characteristic function 
of the sphere: 



2 

(38) F{q) 



d'^r e'^-^ 

r<R 



One can show that F is an integrable function peaked in q = and with a width ~ 1/R. 
From the Parseval-Plancherel identity the integral over the whole momentum space of 
F is {27r)'^Vd{R)- When RSko ^ 1, where Skg is the momentum width of the finite 
temperature Fermi surface, see Eq. ([28|) . we can replace F by a Dirac distribution: 



(39) 



Fig-) {27: fVd(,R)S\q). 
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This leads to 
(40) 



R 



Vd{R) 



-^nkil-Uk). 



To easily calculate the resulting integral to leading order in T, we note that 
(41) d^iik^ (ink{l~nk) 

so that the T = value of the density only is required, Eg. (fT8|) . We finally obtain 



(42) 



2n 



The above result does not apply when RSko ^1. In the range of radius R <^ 
[Sko]^^ we may simply perform a T = calculation. The technique to take advantage 
of the translational and rotational invariance of 0^(|r — r '|) is to introduce the purely 
geometrical quantity 

(43) K[X)= f d'^r [ d'^r'S{\f-f'\- X), 

Jr<R Jr'<R 

SO that 
(44) 



YaTNR = pVdiR) 



2R 



dX K{X)f{X). 



An exact calculation of K{X) is possible: for X < 2R, one obtains 

(45) KiDiX)=A{R-X/2) 

(46) K2DiX)^TrX 



X 



4i?2arccos —-X (4R^ - X^^^ 
2R ^ ' 



(47) 



Kzd{X) = y^'(^ + 4i?)(X - 2Rf. 



When combined with Eas. (|25l27p . this leads to exact expressions for the variance in ID 
and 3D, in terms of the complex integral Ci and Si functions. The large k^R expressions 
are then readily obtained. The asymptotic expansion of the 2D has to be worked out 
by hand, singling out the contribution of the low X quadratic expansion of K{yC). We 
obtain: 

(48) d=l: VariV^ = I±l±|ii^ + 0(l/(fcoi?)^) 

(49) d = 2: YarNii^^ [ln(4fcoi?) + C]+ 0{\n{koR) /k^R) 



(50) d = 3: YavNR^^\n{A^kaR)[koRf-^ 



HBskoR) + 0{l/koR) 
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with 

(51) lim 2^3(^4' 2. i 3; -ii2)_lni.~ 0.656657 

(52) A-i = exp(7 + 21n2- 1/2) = 4.321100... 

(53) ^3 = cxp(7 + 21n2-5/12) = 4.696621... 

and 2^3 is a hypergeometric function, 7 = 0.5772156649 ... is Euler's constant. 

1'4. Application to the ID gas of impenetrable bosons. - Consider a ID homogeneous 
gas of bosons interacting via a deha potential gS{xi — X2) in the hmit g — > +00. In 
this so-caUed impenetrable boson limit, the interaction potential can be replaced by 
the contact condition that the many-body wavefunction vanishes when two bosons are 
in the same point. On the fundamental domain of ordered positions of the N bosons, 
xi < . . . < xn, one then realizes that the eigenwavefunction for bosons coincides with 
an eigenwavefunction of TV non-interacting fermions, with the same eigenenergy. Out of 
the fundamental domain, the bosonic and fermionic wavefunctions may differ by a sign. 

As a consequence, the spatial distribution of the impenetrable bosons, being sensitive 
to the modulus squared of the wavefunction, will be the same as for the ideal Fermi 
gas. The discussion on the variance of the number of particles in an interval [—R, R] in 
subsection ll'3l immediately applies. 

On the contrary, the first order coherence function for the impenetrable bosons (x) 
is sensitive to the phase of the wavefunction and will differ from the gi fermionic function. 
More details, obtained with the Jordan- Wigner transformation, are given in [T] and 
references therein. At zero temperature, we simply recall the mathematical fact, showing 
the absence of condensate, that at large distances [5] 

where kp ~ np, p being the density of bosons, and A ~ 0.92418 . . .. We also use the fact 
(that one can check numerically [T]) that the power law behavior of gi{x) is the same 
as the one of the simpler function 

(55) G{x) ^ (e^'^^i".-!) 

where -/V[o,z^] is the operator giving the number of fermions in the interval from to a; 
(here a; > 0) and the expectation value is taken in the ground state of the ideal Fermi 
gas. In this way, one relates the long range behavior of g^ to the counting statistics of 
an ideal Fermi gas: G{x) is the difference between the probabilities of having an even 
number and an odd number of fermions in the interval [0,a;]. One may assume that 
the probability distribution of the number n of fermions in the interval has a Gaussian 
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envelope [5]. Using the Poisson formula 



+00 



+00 



(56) 



n— — 00 



n— — 00 



where f{k) = J dx f{x) ey:j){~ikx) is the Fourier transform of the arbitrary function /, 
and restricting to the leading terms n = and n ~ 1^ we obtain 



(57) 



G{x) ~ 2cos(7rpa;)e~'^' (Var W[„,.])/2^ 



From the asymptotic expression Eg. ([48)) we obtain 



(58) 



cos(7rpx) 



which explains the decay of gi{x). 



1'5. In a harmonic trap. - In a harmonic trap, simple approximate results can be 
obtained in the limit where the chemical potential is much larger than the trapping fre- 
quencies LUa times 7i, in a temperature range to be specified, in the so-called semi-classical 
approximation. Wc illustrate this semi-classical approximation for two quantities, a ther- 
modynamic one, the entropy, and a local observable, the density. 

1'5.1. Semi-classical calculation of the entropy. The entropy S* is a physically ap- 
pealing way to evaluate to which extent a system is cold, as it is a constant in thcrmo- 
dynamically adiabatic transformations, contrarily to the temperature. 

To calculate thermodynamic quantities in the grand-canonical ensemble, it is conve- 
nient to start from the grand-potential 



(59) 



f}(M,r) 



-/cfiTlnTr e 



-/3(ff-MiV) 



Then one uses the fact that 



(60) 



-SdT-Ndfi 



where N = (N) is the mean number of particles. Since the various field eigenmodes 
decouple in the grand-canonical ensemble, one finds 



(61) 



n{fi, T) = UbT J dep{e) ln[l - n(e)] 



where n(e) is the Fermi-Dirac formula and p{e) is the density of states for a particle in 
the trap: 



(62) 



e - E("a + l/2)?ia;„ 
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where the sum is over all vectors with non-negative integer coordinates. In the semi- 
classical limit one replaces the sum by an integral 



(63) 



J[0.+oo['^ 



where we dropped the ground mode energy eo under the assumption e eo- By the 
change of variables = Uae/huja one obtains 



(64) 



Psc{e) 



-'[0,+oo[<i \ a / 



The integral can be calculated by integrating over Ud and by performing the change of 
variables ri = ui, T2 = ui + U2, ■ • ■ , tj-i = ui + . . . + Ud-i- Wc finally obtain 



(65) 



Psc(e) 



{d-l)\{hujY 



where uj = iX\a ^a)^^'^ the geometric mean of the trap frequencies. 

We replace p by psc in Eq. (|6ip . and we integrate by parts, taking the derivative of 
the ln[l — n(e)], to obtain the semi-classical approximation: 



(66) 



r2 ~ 



de — n{e). 



In the degenerate regime, we use the technique sketched around Eg. pi)) . based on the 
fact that n(e) — nT=o(e) is a narrow function of the energy, to obtain the low-T expansion: 



(67) 



d+1 



{d + iy.ihujy 



1+ d{d + i) 

6 



Calculating 5' = —Dt^ and N = —d^il, we obtain the semi-classical approximation for 
the number of particles and for the entropy per particle: 



(68) 
(69) 

with the Fermi energy 
(70) 



N 



p" 



d\{nujY 



dir'^ 

S/N ~ —kBT/Tp 



keTp = huj{Nd\) 



What are the validity conditions of this approximate formula for the entropy? The 
temperature should be T ^ Tp but, contrarily to the thermodynamic limit case, the 
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temperature should remain high enough, to allow the approximation of the density of 
states (which is a sum of Dirac peaks) by a smooth function. The typical distance 
between the Dirac peaks at the Fermi energy should be smaller than the thermal energy 
width ~ of the Fermi surface. In the case of an isotropic harmonic trap, cua = lo for 
all d dimensions, the Dirac peaks are regularly spaced by fiuj so that we impose 

(71) ksT > huj. 

In the case of incommensurable trap frequencies, the peaks are irregularly spaced and 
we estimate a typical spacing from the inverse of the smoothed density of states at the 
Fermi surface: we require 

(72) keT > ^ 



Psc(/i)' 

with /i ~ ksTp. This leads to the condition 

(73) ,,r»if#.«..(^''"^ 



dN d7Vi-i/d ■ 

One sees that this is a much weaker condition than in the isotropic case, for d > 1, since 
we are here in the large N limit. This fact is illustrated by a numerical example in 3D, 
calculating the entropy for an isotropic trap and an anisotropic trap, see Fig. [21 

1'5.2. Semi-classical calculation of the density. The semi-classical approximation for 
the density profile in the external potential U{r) is strictly equivalent to the so-called 
local density approximation (LDA), as can be revealed by the writings: 

(74) Psc(r) = j J^^Pir) +pV2m] 

(75) Plda(?") = Phom[M - Uif),T] 

where phom is defined in Eq.ljTS]) and n(e) is the Fermi-Dirac formula Eq. pTj) . The 
philosophies of the two approximations of course differ. In the semi-classical case, one 
relies on a phase space density of the particles, which is an approximation to the Wigner 
representation W{f,p) of the one-body density operator of the gas [B]: 

(76) p ^ 



exp[f3{ho - + l 
where the one-body Hamiltonian ho is defined in Eq.(I3]). One assumes 



{2^nY 2 2 ' {2TTnY exp{/3[[/(f ) + 1^ _ + 1 
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0.2 0.4 0.6 0.8 



rJ 

Fig. 2. - Entropy 5 of a spin-polarized ideal Fermi gas in a 3D harmonic trap. Black disks: exact 
result for trap frequencies lj^ = <:i;(0. 5981188 . . . , 1.2252355 . . . , 1.36440128 . . .). Empty circles: 
exact result for an isotropic trap = t^(l, 1, 1)- Solid line: approximate formula Eq.(j69}. The 
chemical potential is fixed to = eo + 47.6?ta), where eo is the trap ground mode energy. In 
order to test the accuracy of Ea. H69p . in a situation where the number of particles rather 
than the chemical potential jj. is known, as it is the case in experiments with atomic gases, we 
use the exact, numerical value of A'^ to evaluate Eq. (|69p rather than its approximation Ea. (|68|l . 
Here iV is about 2 x 10*. 



It remains to integrate this approximation of W{r,p) over p to obtain psc{r). 

In the local density approximation, the gas is considered as a collection of quasi- 
macroscopic pieces that have the same properties as the homogeneous gas of temperature 
T and chemical potential ^lda = I-^ — U{r). Note that this last equation intuitively 
expresses the fact that the chemical potential is uniform in a gas at thermal equilibrium. 

These approximations assume that the external potential varies weakly over the corre- 
lation length of the homogeneous gas, ~ l/ko where fco(^) is the local Fermi wavevector, 
such that h'^kQ/2m = /iLDA- In a harmonic trap, the scale of variation of the potential 
may be taken as the spatial extension of the gas, the so-called Thomas-Fermi length Ra 
along direction a such that 

(78) ^i^^mLolRl. 

The condition ko(f)Ra ^ 1, and correspondingly the semi-classical approximation and 
the LDA, will fail close to the borders of the cloud where fco diverges. 
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At r = 0, as is apparent from Eq. p8[) . these approximations predict the density 
profiles: 



(79) p{f) = 



(2nf 



Integrating this density profile over the domain |rQ,| < Ra gives an expression of the total 
number of particles as a function of the chemical potential: this expression coincides with 
the righthand side of Eg. ([55)1 . whatever the dimension of space d, which illustrates the 
consistency of the various approaches. 

It is instructive to see if the LDA for the density profile can be used at T = strictly, 
or if a lower bound on T is required, as was the case for the entropy. An answer to this 
question is obtained by calculating exactly the density of the gas and its second order 
derivatives in the trap center and comparing to the LDA. The value of p(0 ) is simple 
to obtain from a numerical summation over the harmonic oscillator quantum numbers 
Ua, since the value in a; = of the ID harmonic oscillator eigenstate wavefunction of 
quantum number n > is known exactly: for n odd, (/)n(0) = 0, and for n even, 

where the harmonic oscillator length is a^o = {h/mujY^^. This may be obtained from 
properties of the Hermite polynomials _ff„, using (j)nix) = A/'„-ff„(x/aho) exp(— 2:^/20^^) 
with a positive normalisation factor such that 

-2 = aho7ri/22"n!, or from the recursion 
relation (/)„(0) = —(1 — l/n)^/^0„_2(O), obtained from the a and formalism. The 
second order derivatives of the density in r = are also easily obtained: they involve 
= V^2ri0„-_i(O)/aho and the second order derivative of which is related to (j>n 
by Schrodinger's equation, 0''(O) = — (2n+ !)(/)„ (0)/a^Q. 

We apply this check in ID. At zero temperature, explicit expressions can be obtained. 
For a chemical potential n + 1/2 < fi/huj < n + 3/2, where n is an even non-negative 
integer, we obtain 

2,n^_ 1 



(82) p"(0) = -2a^^p{Q). 

For a chemical potential n + 3/2 < fi/huj < n + 5/2, where n is still even, the a; = 
density assumes the same value (since (f>n+i{0) ~ 0) whereas the second order derivative 
changes sign: 

(83) p(0)^ ' + 



7ri/2ai^^ 2"[(n/2)!]2 
(84) p"(0) = 2a~^pi0). 
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The T = exact value of the density as a function of the chemical potential is 
compared to the LDA in Fig. [3^: the LDA does not reproduce the steps but nicely 
interpolates between them. By using Stirling's formula for the exact result, one finds 
that the relative error in the LDA is 0{1/N), where iV is the particle number. If one 
cleverly approximates the value of a step in p(0) by taking the median value of the 
chemical potential, ^ ~ {n + 3/2)hu!, one even finds that the relative error in the LDA 
vanishes as 0(1/A^^). The LDA, on the contrary, gives a totally incorrect prediction for 
p"{0) (not shown in the figure): one has 

(85) pV(0) = ^ 

7r^<oPLDA(0) 

which tends to zero in the large N limit, whereas the exact result oscillates with a 
diverging amplitude. All this can be understood from the fact that the exact density 
profile at T = is an oscillating function of x, oscillating at smaller and smaller scales in 
the large N limit. This fact can be revealed numerically [7], but more elegantly results 
from the exact summation formula : 

m 

(86) ^(l)l{x) = (m + l)(/)^(x) - yTO(r7r+T)0m+i(a;)0m_i(x). 

fc=0 

The LDA of course misses these oscillations. 

What happens at finite temperature ? We expect that, at fc^T ~ ?itj, these oscillations 
are sufficiently reduced to make the LDA more accurate, while the temperature is low 
enough to allow the use of the T = limit of the LDA. This is confirmed by the ksT = huj 
data in Fig. [31 the plateaus in p{0) as a function of fj. are hardly visible (see a) , so are 
the oscillations of p{x) as a function of cc [3], and the oscillations of p"(0) are so strongly 
reduced that the LDA approximation becomes acceptable (see b). 



2. Two-body aspects of the interaction potential 

In real gases, there arc interactions among the particles. In fermionic atomic gases it 
is even possible to reach the maximally interacting regime allowed by quantum mechanics 
in a gas, the so-called unitary regime, without altering the stability and lifetime of the 
gas. This section reviews possible models that can be used to represent the interaction 
potential and recalls basic facts of two-body scattering theory. We restrict to a three- 
dimensional two spin-state, single species Fermi gas; the p- wave interactions among atoms 
of a common spin state are neglected, usually an excellent approximation away from p- 
wave Feshbach resonances. The s-wave interactions arc on the contrary supposed to be 
strong, the scattering length being much larger than the potential range, in the vicinity 
of a Feshbach resonance. 
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Fig. 3. - Central density of a ID harmonically trapped ideal Fermi gas, in the grand canonical 
ensemble, as a function of the chemical potential fi. (a) Density in the trap center a; = 0: exact 
result for T = (black solid line with steps), exact numerical sum for ksT = Tiuj (red smooth 
solid line), vs the T = LDA (dashed line), (b) Second order derivative of the density in the 
trap center, <f p/dx^{x = 0), for fc^T = hui: exact numerical sum (solid line) vs the T = LDA 
(dashed line). 



2'1. Which model for the interaction potential ? . - The detailed description of the 
interaction of two atoms is involved. At large enough interatomic distances, in particular 
much larger than the size of the electronic orbitals of an atom, one can hope to represent 
this interaction by a position dependent potential V(f)^ which includes the van der Waals 
interaction term 

(87) V{r) ^ 

a simple formula that actually neglects retardation effects and long range dipole-dipole 
magnetic interactions. Forgetting about these complications, we can use Schrodinger's 
equation and the Cq coefficient to construct a length 6, called the van der Waals length, 
that we shall consider as the 'true' range of the interaction potential: 




For alkali atoms, h is in the nanometer range. 

At short interatomic distances, however, this simple picture of a scalar interaction 
potential has to be abandoned, and one has to include the various Born-Oppenheimer 
potentials curves coming from the QED Hamiltonian for two atomic nuclei at fixed dis- 
tance, including all the electronic and nuclear spin states, and the motional couplings 
among the Born-Oppenheimer curves due to the finite mass of an atom. 

Fortunately we are dealing with gaseous systems: the mean interparticle distance is 
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much larger than the potential range b: 

(89) p^/^b < 1 

where p is the mean density. Wc shall also consider atomic dimcrs, but these dimers 
shall be very weakly bound, with a size on the order of the scattering length a, which 
is much larger than b in the vicinity of a Feshbach resonance. As a consequence, the 
central postulate for the theory of quantum gases is that the short-range details of the 
atomic interactions are unimportant, only the low- momentum behavior of the scattering 
amplitude between two atoms is relevant. For fermions of equal masses, this postulate 
has proved robust in its confrontation to experiments, even in the unitary regime. 

A practical consequence of this postulate is that any short range model for the inter- 
action leading to almost the same scattering amplitude fi; as the true interaction, in the 
typical relative momentum range k of atoms in a quantum gas, is an acceptable model. 
We then put the constraint on any acceptable model for the interaction; 

(90) /r'*"' ^ fk 

for the relevant values of the relative momenta k populated in the gas. We insist here 
that we impose similar scattering amplitudes over some momentum range, not just equal 
scattering lengths. Typical values of k can be the following ones: 

(91) fctyp e {a"\fcf,A-i} 

where a is the s-wave scattering length between opposite spin fermions, kp = (Gtt^pi)^/"^ 
is the Fermi momentum expressed in terms of the density of a single spin component 
Pi = Pi, and A is the thermal de Broglie wavelength Eq.([T]). The choice of the correct 
value of ktyp is left to the user and depends on the physical situation. The first choice 
^typ ^ is well suited to the case of a condensate of dimers (a > 0) since it is the 
relative momentum of two atoms forming the dimer. The second choice ktyp ^ kp is well 
suited to a degenerate Fermi gas of atoms (not dimers). The third choice /ctyp ~ A""'^ is 
relevant for a non-degenerate Fermi gas, a case not discussed in this lecture. 

A condition for the gas not to be sensitive to the microscopic details of the potential 
is then that 

(92) A:typ6<l. 

As we shall see, the scattering amplitude in this momentum range should be described 
by a very limited number of parameters, which allows to use very simple models in the 
many-body problem. 

The situation is more subtle for bosons, or for mixtures of fermions with widely 
different masses, where the Efimov phenomenon takes place (see the lecture by Gora 
Shlyapnikov in this volume): at the unitary limit, an effective 3-body attractive inter- 
action occurs and accelerates the atoms; the wavevector k can then, for some type of 
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Feslibach resonances (the broad Feshbach resonances to be defined below), reach the 
range ^ b^^, in which case the details of the true interaction may become important. 
We shall not consider this case here. 

2'2. Reminder of scattering theory. - Suppose for simplicity that two particles of 
mass m interact via the short potential V{ri — that tends to zero at infinity faster 
than l/ri2- Moving to the center of mass frame, with a center of mass with vanishing 
momentum, one gets the Schrodinger equation at energy E: 



(93) E4){r) 



j-2 

-—A,. + v{r 

m 



At negative energies E < 0, one looks for the discrete values of E such that (f> is square 
integrable: these eigenstates correspond to two-body bound states, that we call dimers. 
At positive energies E > 0, one sets 

(94) 

m 

with fco > 0; one looks for scattering states that obey the boundary conditions suited to 
describe a scattering experiment: 

(95) 0(f) = 0o(r) + 0s(r) 

(96) Mr ) = (r |fco) ^ ^^^"'^ 

(97) 0s (f) fkS^^~^ ^ 

The wavefunction 00 represents the free wave coming from infinity, taken here to be a 
plane wave of wavevector fco of modulus fco- The function (/)s(r ) represents the scattered 
wave which, at infinity, depends on the distance as an outgoing spherical wave and pos- 
sibly depends on the direction of observation n = r/r through the scattering amplitude 

It is important to keep in mind that the Schrodinger equation with the above boundary 
conditions is formally solved by 

(98) \4>) = \l + G{E + iQ+)V]\M 

where G{z) — (z — H)^^ is the resolvent of the Hamiltonian, z a complex number. Using 
the identity G = Go + GqVG, where Go{z) = [z — Hq)^^ is the resolvent of the free, 
kinetic energy Hamiltonian, one gets 

(99) 10) = [1 + Go(S + iO+)TiE + iO+)] |0o) 
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where we have introduced the T operator (or T-matrix) 
(100) T{z) = V + VG{z)V. 



Eq. (|99p allows to obtain the equivalence in momentum space of the asymptotic conditions 
Eas. (|95l96l97p in real space: 

(101) {k\^) = {27rfSHk-ko) + ^ , ^ ^2,2/ + ^0+)\ko), 

where the matrix element of T is an a priori unknown smooth, regular function of k. The 
link between the /c-space and the r-space points of view is completed by the following 
relation: 

(102) f.^(^ri)^-^{kon\T{E + iO+)\ko), 

which follows from the expression of the kinetic energy Green's function for E > 0: 



(103) {r\G^{E + iQ+)\r') = - 



Clearly, the central object is the scattering amplitude. It obeys the optical theorem 

(104) a,eatt(fco) = J d^n \ff:jn)\^ = ^Im fj:^{k, / ko) , 

where CTscatt is the total scattering cross section for distinguishable particles. This the- 
orem is a direct consequence of the conservation of probability: the matter current J in 
the stationary state 4> has a vanishing divergence, that is a zero flux through a sphere of 
radius r — > +oo: using the asymptotic form of 0, one gets the announced theorem. Here 
we shall consider model potentials that scatter only in the s-wave: (/j^ is isotropic and so 
is the scattering amplitude, which reduces to an energy dependent complex number: 

(105) f^^{n)^h,. 
The optical theorem simplifies to 

(106) l/fco 



2 _ Im fko 



ko 

Using Imz^^ = — Imz/jzp, for a complex number 2, one realizes that this simply im- 
poses: 
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where u(ko) is at this stage an arbitrary real function. The scattering amphtude takes 
its maximal modulus when u is negligible as compared to fco, this constitutes the unitary 
limit 

A last point concerns the two-body bound states: their eigenergies correspond to poles 
of the resolvent G, which are also poles of the T matrix. A simple way to find the dimer 
eigenenergies is therefore to look for poles in the analytic continuation of the scattering 
amplitude to negative energies, setting fcp = iqo, qo > 0, so that E = — ?i^g§/m < 0. Note 
that the determination +iqa is chosen for fco = VmE/h so that the matrix elements of 
Go, see Ea. (|103|) . assume the correct value (that tends to zero at infinity). Furthermore 
one can easily have access to the asymptotic behavior of the dimer wavefunctions (jiji 
including the correct normalization factor. One uses the closure relation 

(109) / ^l^£o)('^4l+El^.-)(^.l=^ 

where / is the identity operator. One takes the matrix elements of this identity between 
(r I and |r'), for large values of r and r' so that the asymptotic expression Eg. iP?)) may 
be used. After angular average (for a s-wave scattering), one gets factors e='='''o (''+'' ^ and 
^iko{r-r ) ^ Using the optical theorem one shows that the coefficient of e*'^o(''~'' ) vanishes. 
Assuming that u{k) is an even analytical function of k one can extend the integral over 
fco to ] — cx),+oo[, setting f^^ = f-kg, and use a contour integral technique by closing 
the contour with an infinite half-circle in the Imz > part of the complex plane. From 
Eq. (|107[) the residue of fk in the pole k = iqj is obtained in terms of u'{iqj). We finally 
obtain the large r behavior of any s-wave correctly normalized dimer wavefunction: 



(110) Mr) 



1/2 -or 

qj \ ' p. 'ii^ 



l-iu'{iqj)J y/2Kf 



2'3. Effective range expansion and various physical regimes. - For the quantum gases, 
it is very convenient to introduce the usual low-fc expansion of the scattering amplitude: 

^^^^^ " ' a-^+iko-klrj2 + ... 

where the parameter a is called the scattering length, the parameter Tg is called the 
effective range of the interaction, not to be confused with the true range b. As discussed 
below, this expansion is interesting if the omitted terms ... in the denominator of Eg. ljllip 
are indeed negligible when Eg. ([5^ is satisfied. 

We take the opportunity to point out that the general allowed values for the effective 
range can go from — oo to -l-oo. This can be demonstrated on a specific example, the 
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square well potential, with V{r) = —Vq for r < b, V{r) = otherwise, the true range 
b being fixed and Vq being a variable positive quantity. Setting kq = (mVo/h^)^^^ , we 
obtain the exact formulas: 

(fl2) a^b-'-^ 

Kq 

This easily shows that for all non-zero values of kq such that the scattering length van- 
ishes, the effective range tends to — oo. When kq tends to zero, the scattering length 
tends to zero (more precisely, 0"), quadratically in Kg, as is evident also from the Born 
approximation; in the expression for r^, the last two terms diverge to infinity, with op- 
posite sign; a systematic series expansion in powers of Kob leads to 

(114) re = -— + -^& + 0(Kg&3) forKo^O+. 

5a 175 

One then sees that tends to -l-oo for kq 0"*" since a ^ 0~ in this limit. Note that 
the leading term in Eq. (|114p can also be obtained in the Born approximation. 

As already mentioned, a very favorable case is when the ... in the denominator 
of Eq. (|llip is negligible as compared to the sum of the first three terms, in the low 
momentum regime fco^ ^ 1- one can then characterize the true interaction by only two 
parameters, the scattering length and the effective range. In what follows, it is implicitly 
assumed that this favorable case is obtained, otherwise more realistic modeling of the 
interaction should be performed. One can then identify interesting limiting cases. 
The zero-range limit: This is a limit where the only relevant parameter of the true 
interaction is the scattering length a, i.e. we shall assume that the typical momentum 
ktyp satisfies 

(115) fc2yp|re| < |a"i -fifctypl, 

where r^ is the effective range of the true interaction. In practice, in present experiments 
and for typical atomic densities, this is true for the so-called broad Feshbach resonances, 
where r^ ~ b: the condition Eq. (|115p is then a direct consequence of Eg. ([5^ . For narrow 
Feshbach resonances, |re| can be considerably larger than b and Eq. (|115p may be violated 
for typical densities. 

What type of model can be chosen in this zero-range limit ? A natural idea is to 
use a strictly zero-range model, having both a vanishing true range and effective range, 
and characterized by the scattering length as the unique parameter. This 'ideal world' 
idea corresponds to the Bcthe-Peierls model of subsection 12 51 where the interaction is 
replaced by contact conditions. 

In practice, the contact conditions are not always convenient to implement, in an 
approximate or numerical calculation. It is therefore also useful to introduce a finite 
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range model, usually with a true range of the order of its effective range r™° . One 
then usually puts the model range to zero at the end of the calculation. Note that this 
is not always strictly speaking necessary: it may be sufficient that the model is also in 
the zero-range limit, 

(116) kl^\rrn<^\a'' + ^hy^\- 

In section [3] on BCS theory we shall restrict to the zero-range limit. We shall then 
approximate the true interaction by an on-site interaction in a Hubbard lattice model, 
see subsection 12 61 the effective range of the model is then of the order of the lattice 
spacing, and one may easily take the limit of a vanishing lattice spacing in the BCS 
theory. 

The unitary limit: What is the condition to reach the unitary limit in a gas, in a regime 
where the low-A: expansion holds ? One should have 

(117) |a-i|,/c2yp|re|«fctyp. 

The condition on a is in principle possible to fulfill, as a diverges close to a Feshbach 
resonance. The condition on Tg shows that the zero-range limit Eq. piSp is a necessary 
condition to reach the unitary limit. This is why the broad Feshbach resonances are 
favored in present experiments (over the narrow ones) to reach the strongly interacting 
regime. To make it simple, one may say that the unitary limit is simply the zero-range 
limit with infinite scattering length. 

2'4. A two-channel model. - We have the oversimplified view of a magnetically induced 
Feshbach resonance, see Fig. ID two atoms interact via two potential curves, Vopen{r) and 
V^iosod(?')- At short distances these two curves correspond to the electronic spin singlet 
and triplet, respectively, for atoms of electronic spin 1/2, so that they experience different 
shifts in presence of an external magnetic field. At large distances, Vopen('') tends to zero 
whereas Vciosod(f) tends to a strictly positive value Voo- 

The atoms are supposed to come from infinity in the internal state corresponding to 
the curve Vopen(''), the so-called open channel. The atoms are ultracold and have a low 
kinetic energy, 

(118) E<^Voo. 

In real life, the two curves are actually weakly coupled. Due to this coupling, the atoms 
can have access to the internal state corresponding to the curve T4ioscd(?'), but only 
at short distances; at long distances, the external atomic wavefunction in this so-called 
closed channel will be an evanescent wave that decays exponentially with r since E < Voo ■ 
Now assume that, in the absence of coupling between the channels, the closed channel 
supports a bound state of energy -Emoi, denoted in this text as the molecular state. 
Assume also that, by applying a judicious magnetic field, one sets the energy of this 
molecular state close to zero, that is to the dissociation limit of the open channel. In 
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this case one may expect that the scattering amphtude of two atoms may be strongly 
affected, by a resonance effect, when a weak coupling exists between the two channels. 



> 




Fig. 4. - Oversimplified view of a magnetically induced Feshbach resonance. The interaction 
between two atoms is described by two Born-Oppenheimer curves. Solid line: open channel 
potential curve. Dashed line; closed channel potential curve. When one neglects the coupling 
A between the two curves, the closed channel has a molecular state of energy i?moi with respect 
to the dissociation limit of the open channel. Note that the energy spacing between the solid 
curve and the dashed curve was greatly exaggerated, for clarity. 



A quite quantitative discussion is obtained thanks to a very simple two-channel model. 
This model is the most realistic of the models presented here, and we shall consider it 
as a 'reality' against which the other models should be confronted. Assuming that the 
particles are free in the open channel, and can populate only the molecular state in the 
closed channel, the only non-trivial ingredient is a coupling between the molecule with 
center of mass momentum zero, |mol), and a pair of atoms of opposite momenta and spin 
states ||fc), with 

(119) ||fc)^a[,y_^.jO), 

where the creation and annihilation operators obey the free space anticommutation re- 
lation 



(120) 
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This coupling in the momentum representation is defined as 
(121) V\mo\) = [ -^Ax(fc)llfc) 



^ (27r)3 

(122) V\\k) ^ Ax{k)\mo\). 

We have introduced a real coupling constant A, and a real and isotropic cut-ofF function 
x(fc), equal to 1 for fc = and to for k = oo, to avoid ultraviolet divergencies. It is 
convenient to take 



(123) x{k) 



e 



where the length ct is a priori of the order of the true potential range b. The two-channel 
model is thus defined by the three parameters £',„oi,A and a. It scatters in the s-wave 
channel only, since x is isotropic. 

An eigenstate of energy E of the two-channel Hamiltonian is of the form 

j^a{k)\\k). 

Including the molecular state energy and the kinetic energy of a pair of atoms, the 
stationary Schrodingcr equation leads to 

(125) a{k) + Ax{k)P = Ea{k) 

m 



(126) i?n,oi f3 + / 7— 3 Ax(fc)a(A; ) = Ep. 



(2^ 



We restrict to a scattering problem, E > Q. Instructed by the previous discussion around 
Eg. pOip . we correctly solve for a in terms of /3 in Ea. (|125p : 

(127) a{k) = {2^)H{k ko) + ^ , .^t^^^l,,, /3, 

E + iO+ — n k^/m 

where ko is the wavevector of the incident wave, and E — h^kQ/ra. Insertion of this 
solution in Eq. (|126p leads to a closed equation for (3 that is readily solved. From Eq. (jlOip 
and Ea. (|102p we obtain the scattering amplitude 

(128) - 



'inn. - + J j2Fj^ E+tO+-h-^kym 



Let us work a little bit on this scattering amplitude. First, using the identity 
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where V denotes the prmcipal value, one can check that fk^ is indeed of the form Eq. (|107p . 
so that it obeys the optical theorem. Second, the scattering length is readily obtained 
by taking the fco ^ limit: 



(130) 



1^ 



(Pk mx^ik) 



(27r)3 



This shows that the location of the resonance (where = 0) is not i?moi = but a 
positive value -E^joi, this shift being due to the coupling between the two channels. For 
^-moi larger than this value, a < 0. For -Emoi below this value, a > 0. 

Another general result, not specific to a Gaussian choice for x(fc) (but assuming that 
X(fc) varies quadratically around A: = 0), concerns the zero range limit ct ^ for fixed 
values of fco, A and the scattering length a (note that the molecular energy is then not 
fixed): 



(131) 



/ko 



+ ifco + 



E 



This can be seen formally as resulting from the identity °° dX VI / {X^ — 1) = 0. This 
shows that the two-channel model has a well defined limit, in the zero true-range limit, 
characterized only by the scattering length and by a non-zero and negative effective 
range, 



(132) 



This limit may be described by generalized Bethe-Peierls contact conditions [lOJ. This 
illustrates dramatically the difference between the true range and the effective range of 
an interaction. This also shows that, in this limit, the two-channel model can support 
a weakly bound state, that is a dimcr with a size much larger than the true range a. 
Setting fco = igo, one finds that the right hand side of Ea. p3ip has a pole with go real 
and positive if and only if a > 0: 



(133) 



90 



1 - Vl - 2r0/c 



Finally, an exact calculation of the scattering amplitude for the Gaussian choice 
Eq. ()123p is possible: the trick put forward by Mattia Jona-Lasinio is to calculate the 
analytic continuation to imaginary fco = igo, so that one no longer needs to introduce 
the principal part distribution. One gets 



(134) 



(/igo ) 



qoeric{qoa), 
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where erfc is the complementary error function, that tends to unity in zero, erfc(x) 
1 — 27r~^/^a; + 0{x^). The scattering length is given by 

The calculation of the effective range for a finite a is then straightforward: 
(136) re=rl + ^(l- 



2a 

We shall then classify the Feshbach resonances as follows: for a = oo, the ones with 
Te ~ < 0, much larger than the true range in absolute value, constitute narrow 
Feshbach resonances. The ones with — a constitute broad Feshbach resonances. 

This classification is not only theoretical but also reflects the magnetic width AB of 
the resonance. If one assumes that i^moi is an affine function of the magnetic field B, 
with a slope — /x, one finds from Ea. (|130p that the scattering length a of the two-channel 
model varies with B as 

where Bq is the magnetic field value of the resonance center, and we have singled out 
the value Obg of the background scattering length (that is the value of the scattering 
length of the true interaction far from the resonance), motivated by the fact that a = 
flbg[l + A_B/(i?o — B)] in a more complete theory including the direct interaction in the 
open channel [TT]. This defines the width AB, such that 

(138) ^lAB= ^. 

So, when Obg ^ a ^ b, comparing ^AB to the natural energy scale /mb^ amounts to 
comparing jr^j to b. 

2'5. The Bethe-Peierls model. - In this model, the true interaction is replaced by 
contact conditions on the wavefunction. For two particles in free space, in the center of 
mass frame, the contact condition is that there exists a constant A such that 

(139) <j3{r) = A[r-^ - a-^] +0{r) 

in the limit r 0. For r > the wavefunction is an eigenstate of the free Hamiltonian: 

(140) E(j)(f) = Ap(j)(f). 

m 
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In mathematical terms, this defines the domain of the Hamiltonian. The Hilbert space 
remains however the same, with the usual scalar product, because a 1/r divergence in 
3D is square intcgrablc. 

An equation valid for all values of r can be written using the theory of distributions: 

. , 'iiTh^a 

141 E^{r) = A,-0 f + cbres S{r) 

m m 

where the so-called regular part of (jj is 

(142) 0„g = - dr N(f )],,^o = -a-^ lim r0(f ). 

a r— tO 

This establishes the equivalence with the regularized pseudo-potential model [TH [HI HH 
I15j . The last form of (preg is convenient to use in the limit a = oo, since simplifies 
with the factor a in front of the Dirac S{f). 

The solution of Ea. p40p with the boundary conditions Ea. (|W|) and the contact con- 
dition is simply: 

(143) ^^r^) ^ e^ka-r + fk/-— 

r 

with the scattering amplitude 

(144) fk„ 



a ^ + ik 

This shows that the model has both a zero true range and zero effective range. It may be 
applied to mimic the effect of the true interaction potential when the condition Eq. (|115p 
is satisfied. 

Following the discussion above Eg. ljllOp . one finds that the model supports a dimer 
if and only if a > 0. The pole of fiq^ is then go = l/^j resulting in a dimer energy 



(145) Eo = - 



Since the model has a zero true range, it turns out that the dimer wavefunction coincides 
everywhere with its asymptotic form Eq. (jllOp , 



1 

(146) Mr) 



'2'Ka r 

The advantage of the Bcthe-Peierls model for the A^-body problem is that it introduces 
the minimal number of parameters (the scattering length). In the unitary limit, it is a 
zero-parameter model, which is ideal to describe universal states. This advantage actually 
really exists if one restricts to analytical solutions of the model (numerical solutions 
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tend to introduce extra parameters like a discretization step or a momentum cut-ofF). 
Analytical solutions in free space are available up to = 3 in the unitary limit [HI El ■ 
For particles trapped in harmonic potentials, analytical solutions arc available for two 
particles at any a for isotropic traps [16j and cylindrically symmetric traps |17| . For 
three particles in an isotropic harmonic trap, the problem was solved analytically only 
in the unitary limit a = oo [501 HI] ■ In the many-body problem, for the universal unitary 
gas in an isotropic harmonic trap, a scaling solution can be found for the many-body 
wavefunction in the time dependent case |22| and an exact mapping to free space zero 
energy cigcnstates can be constructed [531 (211 US . 

It is worth mentioning that, for N — i bosons, the Bethe-Peirls model docs not define 
a self-adjoint operator and has to be supplemented by extra contact conditions; this is 
related to the Efimov phenomenon [THl[Tn]; the problem of course persists for TV > 3, 
making the model non satisfactory. Fortunately this lecture is devoted to equal mass 
spin 1/2 fermions, for which no Efimov effect appears; the corresponding unitary gas is 
then called universal, to emphasize this aspect that no extra parameter has to be added 
to the Bethe-Peirls model. Note that this absence of Efimov effect is not simply due to 
the fact that the Pauli principle prevents one from having three fermions of spin 1/2 in 
the same point of space: if the spin up and spin down fermions have widely different 
masses, an Efimov effect may appear |26| . 

We take the opportunity to mention another trap to avoid, even for fermions. With 
the formulation of the model in terms of binary interactions with the pseudo-potential 
Vpp{_r) = g5{f)dr{r •), where g = 4Trh^a/m, it is easy to 'forget' that contact conditions 
are imposed on the many-body wavefunction. One may then be tempted to perform 
a variational calculation with a iV-body trial wavefunction which docs not satisfy the 
contact conditions Eq. (|139p . that is which is not in the domain of the Hamiltonian [7f| . 
Whereas in the low a-limit the corresponding prediction for the energy may make sense 
as a perturbative prediction, it may become totally wrong in the large a limit. 

To illustrate this statement, let us consider two identical particles in a harmonic 
isotropic trap interacting via the pseudo-potential. After separation of the center of 
mass coordinates, one is left with the following Hamiltonian for the relative motion 



where fj, = m/2 is the reduced mass. Let us take as a trial wavefunction the Gaussian 
of the ground state of the harmonic oscillator. 



(147) 




(148) 



Mr) 



e-'-V(4aL) 



(27r)3/44; 



where Oho = y^h/muj. This trial wavefunction docs not obey the contact conditions, so 
it is not in the domain of the Hamiltonian. It is therefore mathematically incorrect to 
calculate the mean energy of (jjt by representing the action of H on by the differential 
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operator Eq. (|147[) . One rather has to calculate the overlap {ipn\4't) of \(j)t) with each 
exact eigenstate i/'„ of H of eigenenergy En, and calculate the sum J2n ^n\{'ipn\(t>t)\'^ ■ 
From [T6| and for a finite a one finds that the overlap squared scales as 1/ n'^/^ whereas 
the energy i?„ scales as n, for large n. The correct mean energy of (j)t is thus infinite. 
If one however falls in the trap, one gets the wrong mean energy 

(149) E^,ong = (0t|i^|0t) = + ( z ) 



2 \ TT y flho 

In the low a limit this can be compared to the expansion of the exact minimal positive 
eigenenergy (16j 



3 

(150) -Ecxact = + 



0\ 1/2 o / 



+ -(1 -ln2 

TT / flho TT ' \ flho 



One sees that the wrong variational calculation gives the correct linear correction in a, 
that is the correct perturbative result. One may then be tempted to use the wrong 
variational calculation to predict the sign of the coefficient of a^; one would naively 
assume -Eexact ^ ^'wrong and one would wrongly predict a negative sign for the coefficient 
of a? |2H] • In the opposite limit of an infinite scattering length, one directly realizes the 
absurdity of the assumption inexact < -E'wrong: for a — oo, one finds E'wrong ^ — oo, 
whereas -Eexact Tiuj/2, with the exact wavefunction oc e^*" jr. 

Note that the same pathology occurs for the many-body problem. Consider for exam- 
ple a two spin-component spatially homogeneous Fermi gas, with interactions modeled 
by the pseudo-potential. If one uses the Hartree-Fock variational ansatz, which is not in 
the domain of the Hamiltonian, one correctly obtains, in the weakly attractive regime 
kpa — > 0~, the mean field shift of the ground state energy, but one wrongly predicts a 
negative sign for the coefficient of in the low kpa expansion of the ground state energy, 
by wrongly arguing that the exact ground state energy has to be below the 'variational' 
Hartree-Fock energy. As can be checked on the systematic expansion the coefficient 
of is actually positive, because 11 > 2 In 2. 

One should not leave this section with the impression that any variational calculation 
is impossible within the Bethe-Peierls model. To be safe, one simply has to take a 
variational ansatz which is in the domain of the Hamiltonian. An example of a successful 
variational calculation is the derivation of the virial theorem for a unitary gas in a non 
necessarily isotropic harmonic trap; this derivation, due to Frederic Chevy, is reported 
in |25| . and uses the fact, for the universal unitary gas 1/a = 0, that the function of the 
rescaled coordinates \I'(ri/A, . . . , rjv/A), where A > 0, is in the domain of the Hamiltonian 
if the wavcfimction ^'(ri, . . . , rjv) is in the domain of the Hamiltonian. 

2'6. The lattice model. - The last model we consider is at the intersection of two very 
popular classes of models in condensed matter physics, the separable potentials and the 
lattice models. It is slightly simpler than the two-channel model, but it does not apply 
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to narrow Feshbach resonances in a situation where the effective range term cannot be 
neglected in the scattering amplitude. It is more tractable than the Bethe-Peierls model 
for variational treatments (like the BCS method, see section [3|) or for exact numerical 
treatments since the domain of the Hamiltonian is the one of usual quantum mechanics. 
In particular, since it is a sort of Hubbard model, the machinery of quantum Monte Carlo 
methods applicable to the Hubbard model [SHI [SD 1311 could be apphcd to it. 

In this model, the spatial coordinates r of the particles are discretized on a cubic 
grid of step Z. As a consequence, the components of the wavevector of a particle have a 
meaning modulo 27r/Z only, since the hmction f —^ exp^ik ■ r) defined on the grid is not 
changed if a component of k is shifted by an integer multiple of 27r/L We shall therefore 
restrict the wavevectors to the first Brillouin zone: 



(151) kev^ . 

This also shows that the lattice structure in real space provides a cut-off in momentum 
space. 

The interaction between opposite spin particles takes place when two particles are at 
the same lattice site, as in the Hubbard model. In first quantized form, it is represented 
by a discrete delta potential: 

(152) ^=f<5,-,A- 

The coupling constant go is a function of the grid spacing I. It is adjusted to reproduce 
the scattering length of the true interaction. The scattering amplitude of two atoms on 
the lattice with vanishing total momentum is given in [34j , and a detailed discussion is 
presented in [T]. We give here only the result, generalized to an arbitrary even dispersion 
relation k E-^ for the kinetic energy on the lattice: 

(153) /.o = -T^ ^ 



Adjusting go to recover the correct scattering length then gives 

1 1 f d?k I 

(154) 



50 g U{27r)3 2Ej:' 

with g = ATrfi^a/ra. This formula is reminiscent of the technique of rcnormalization of 
the coupling constant. 

A natural case to consider is the one of the usual parabolic dispersion relation, Ej: = 
fi^k^ /2m. A more explicit form of Ea. (|154p is then 



(155) 



AirTi'^a/'m 
~ 1-Ka/l 
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with a numerical constant given by 

r.7r/4 



(156) 



12 f ' 

K ^ — ^6* ln(l + 1/ cos^ 61) = 2.442 749 607 806 335 ... , 

Jo 



and that may be expressed analyticaUy in terms of the dilog special function. The 
effective range of the model is easily calculated with the complexification technique fco = 
iqo; it is positive and proportional to I: 

(157) r^°°d<=i^ / — = 0.336 868 47... L 

27r-^ Jfe^[-l/2,l/2P « 

In the I limit, for a fixed a (not fixed go) the lattice model then reproduces the 
scattering amplitude of the Bethe-Peierls model, 

(158) lim/fc„ =- , 

and admits a dimcr for a > 0. Note that go indeed assumes negative values that tend 
to 0~ linearly wi th I in this limit: the interaction is attractive, whatever the sign of the 
scattering length 

Studying the weakly interacting regime: In the opposite case of / 3> |a|, the bare coupling 
constant is 



(159) 50 



so that go has the sign of the scattering length. For a > 0, the model interaction is 
repulsive, for a < it is attractive. In both cases, it does not support any two-body 
bound state. Its scattering length can be calculated in the Born approximation, since 
the term in a in the denominator of Eq. (|155p is small as compared to unity. This makes it 
an appropriate model to study the weakly interacting regime, as was done for Bose gases 
in [34j : the Hartrec-Fock theory, for example, which relics on the Born approximation, 
may be applied, and its accuracy may be supplemented by keeping higher order terms 
in a perturbative expansion of the interaction. In the context of fermions, it may be 
used in the weakly attractive or weakly repulsive regime; in the latter case, it may be an 
alternative to the celebrated hard sphere model (with hard spheres of diameter a > 0). 

It is instructive to apply the philosophy around Eq. (j90p to understand why the lat- 
tice model with I 3> |a| is restricted to the weakly interacting regime, for a quantum 
degenerate gas of atoms with ktyp — kp. 



(^) More generally, the lattice model admits a dimer if go < g^, where gg is the a ^ oo limit of 
Ea. (|155|l . This condition is equivalent to go < and a > 0. 
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In an effort to have f™°'^°^ ~ fk, one first adjusts the scattering length of the model 
to the correct value a by tuning g^. 

If one further adjusts the lattice spacing I to have r™"*^"' = re- this assumes a 
broad Feshbach resonance with re ~ b, since the narrow Feshbach resonance has 
Te < 0; one then finds / ~ b; since fctyp6 1 (otherwise microscopic details of the 
interaction will come in, and the ... in Eq. ()llip is not negligible), having \a\ ^ I 
leads to fci?|a| <^ 1, which is the weakly interacting regime. 

Alternatively, one may consider the more favorable case of a true interaction being 
in the zero-range limit Eg. pisp . The model should also be in the zero-range limit 
Eq. pi6p , which can be written here 



(160) kpl-t: 



1 



and this is sufficient (r™°'^''' ex I may differ a lot from re). Since we wish to have 
\a\ <IC I, the above equation leads to kp\a\ <§; \{kFa)~^ + i\, a. condition that may 
be satisfied only in the weakly interacting limit /cj^jaj ^ 1. 

• The same result may be obtained in the following much faster way: since we are 
dealing here with fermions, there will be at most one fermion per spin component 
per lattice site, so that kpl < (Btt^)^/'^; the condition \a\ <C I then immediately leads 
to kp\a\ < 1. This reasoning is however specific to the lattice model, whereas the 
previous reasoning is easily generalized to the hard sphere model (where r™°'^''' cx 
a). 

The 'true' Hubbard model: To be complete, we consider a second case for the dispersion 
relation, the one of the true Hubbard model: this makes the link with condensed matter 
physics, and this also shows that a universal quantum gas may be described by an 
attractive Hubbard model in the limit of a vanishing filling factor. The Hubbard model 
is usually presented in terms of the tunneling amplitude between neighboring sites, here 
t = —h'^/2mP, and of the on-site interaction U = ga/l^- The dispersion relation is then 

(161) A^^[l-cos(fc„Z)] 



where the summation is over the three dimensions of space. Close to fc = 0, the Hubbard 
dispersion relation by construction reproduces the free space one h^k^ /2m. The explicit 
version of Eq. (|154p is obtained from Ea. (|155p by replacing the numerical constant K 
by ii-Hub = 3.175911 .... At the unitary limit this leads to U/t = 7.913552 . . ., which 
corresponds to an attractive Hubbard model since t < 0, lending itself to a Quantum 
Monte Carlo analysis for equal spin populations with no sign problem [30l [3TJ [32] . We 
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have also calculated the effective range of the Hubbard model, which remarkably is 
negative: 

(162) rf""^ ~ -0.3056L 

2'7. Application of Bethe-Peierls to a toy model: two macroscopic branches. - In this 
subsection, we move to the problem of N interacting fermions, N/2 with spin component 
t and N/2 with spin component |. In order to have a global view of what the BCS 
theory will be useful for, it is instructive to start with a purely qualitative model |46| . 

Consider a matter wave of isotropic wavcfunction (p{r) in a hard wall spherical cavity 
of radius R, in presence of a point-like scatterer of fixed position in the center of the 
cavity. The effect of the point-like scatterer is characterized by the scattering length a 
and is treated by imposing the Bethe-Pcicrls contact condition Ea. (|139p . 

What is the link between this model and the many-body problem of N interacting 
fermions ? The wavefunction (j){r) describes the relative motion of a, let us say, spin t 
atom, with respect to the nearest spin | atom modelized by the point-like scatterer. The 
cavity represents (i) the interaction effect of the other N/2 — 1 spin J, atoms and (ii) the 
Pauli blocking effect of the other N/2 — 1 spin t atoms. Interaction effect (i): the radius 
of the cavity should then be of the order of the mean interparticle separation in the gas, 

(163) Roc^. 

hp 

Pauli blocking effect (ii): in the case g = 0, the zero point energy of should be of the 
order of the Fermi energy so that one has also the choice (|163p . This explains also the 
choice of hard walls for the cavity; in the case of bosons, imposing that the derivative 
drip vanishes vn r = R [47j . or taking a cubic cavity with periodic boundary conditions 
would be appropriate. Finally, the total energy of the gas is related to e by 

(164) E oc Ne. 

Specific but arbitrary proportionality coefficients are used in Eas. (|163ll64p to produce 
the Fig O to come, as detailed in [15] . 

Let us proceed with the calculation of the eigenenergies of the matter wave in the 
cavity. An cigcnmodc of the cavity with energy e solves Schrodingcr's equation for 
< r < i? 

(165) A0(r) = e</)(r) 

771 

with the Bethe-Peierls contact condition Ea. p39p . If e > 0, we set e = h^k"^ /m and k 
solves 



(166) 



tan kR = ka. 



34 



YVAN CASTIN 



If e < 0, we set e = —h^K^/m and k solves 

(167) tanhfti? = na. 




Fig. 5. - In the toy model with a scatterer in a hard wall spherical cavity: (a) energy per particle 
and (b) pressure of the gas in the first two branches, as functions of —l/kpa. kp, Ep and Pf 
are respectively the Fermi wavevector, the Fermi energy and the pressure of the T = ideal 
Fermi gas with the same density as the interacting gas. 



In the cavity, there is an infinite number of discrete modes. To each of this mode 
we associate a distinct macroscopic state of the gas. In figure [S] we have plotted the 
energy per particle and the pressure P = —dyE, where V = N/p is the total volume of 
the gas, for the first two branches, as functions of —1/kpa. We have taken ~\/kpa as 
the abscissa because it allows an almost direct mapping with the B field axis in a real 
experiment around the location Sq of the Feshbach resonance, see Ea. (|137p . 

The first excited branch is metastable. It starts with a weakly repulsive Fermi gas on 
the extreme left of the figure and has a larger energy than the ideal Fermi gas, indicating 
effective repulsion. When < a <C fc]^^, three-body collisions (not included in the toy 
model) lead in a real gas to the formation of dimers : in the language of the toy model, 
the system starts populating the ground branch. This suggests a first experimental way 
to produce a gas of dimers. 

The ground branch continuously connects the weakly attractive Fermi gas (on the 
right) to a gas of dimers (on the left). This provides a second experimental way to 
produce a condensate of dimers, by adiabatic following of the ground branch. The sharp 
decrease of the total energy on the left part of the figure reflects the 1 /a^ dependence of 
the dimer binding energy eq- The pressure is always less than the Fermi pressure of the 
ideal Fermi gas, indicating effective attraction with respect to the ideal Fermi gas. Note 
that P drops very rapidly on the left side, due to the absence of interaction between 
the molecules in the toy model. In real life, this interaction occurs with a scattering 
length that was calculated by solution of the 4-body problem [48] . As discussed in other 
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lectures of this volume, the molecular condensates have been obtained experimentally 

[asiiioiHiiiia. 

The regime of infinite scattering length \a\ — +00 is universal in the toy model: one 
finds that the energy of the corresponding unitary gas has to be proportional to the 
ground state energy of the ideal Fermi gas: 

(168) ^unitary ^ ^^^idoal 

where E^'^'^^ ~ 3NEf/5, r] is a numerical constant depending on the branch, and Ep 
is the Fermi energy of the ideal gas. In the present state of the art of the field, it is 
accepted that this universality property also holds for a real ground state Fermi gas. 
Approximate fixed node Monte Carlo calculations [351 13S] give for the ground branch 
r] ~ 0.4, in disagreement with early experiments |37| but in good agreement with recent 
experimental measures [3S1 [511 US] and with recent exact Quantum Monte Carlo results 

The toy model allows to easily solve the following paradox: 

• A common saying is that a gas with a positive scattering length a > experiences 
effective repulsion, and that it experiences effective attraction for a < 0. 

• Another saying is that for \a\ ~ +00 the energy of the gas is universal and does 
not depend on the sign of a. 

• However, if a ^ +00 (scenario 1), one expects that the universal state has effective 
repulsion, whereas for a —00 (scenario 2), one expects that the universal state 
experiences attraction. How can it be that there is no dependence on the sign of a 
in the unitary limit? 

The answer provided by the toy model is simple: scenario 1 and scenario 2 arc predicted 
to lead in the unitary limit to two different universal states, belonging respectively to 
the first excited branch and to the ground branch, one experiencing effective repulsion, 
the other effective attraction. 

3. Zero temperature BCS theory: study of the ground branch 

In this section, we consider the many-body problem of a two-component Fermi gas 
with an interaction characterized only by the s-wave scattering length a. 

To be able to use the BCS theory for all values of kpa, while getting sensible results, 
we restrict to the zero temperature case: the usual static BCS theory is indeed unable to 
get a fair approximation to the critical temperature on the bosonic side of the Feshbach 
resonance (a > 0, kpa <C 1), for reasons that will become clear at the end of this section. 
An improved finite temperature theory can be found in [441 145| . 

For simplicity we also assume that the two spin states t and | have the same number 
of particles, so that they have a common chemical potential /i. The case of imbalanced 
chemical potentials and particle numbers is presently the object of experimental and 
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theoretical studies, and may lead to a variety of observed and/or predicted phenomena, as 
the not yet observed non-demixed BCS phases with spatially modulated order parameter 
(the so-called LOFF phases) [49l [50] or the already observed spatial demixing of the two 
spin components [STJ [5ll [S3] . 

Going beyond the toy model of previous section, the zero temperature BCS theory 
predicts that the particles arrange in pairs and gives several properties of the gas of pairs: 

• The existence of long range order in the pair coherence function: 

where ^'(r) = ijj^{r)Tp-^(r) annihilates a pair of particles with opposite spin in 
f. In the thermodynamic limit, the zero temperature BCS theory predicts that 
g^^" has a non-zero limit for large r: the pairs form a condensate. Except in the 
regime kpa — > 0+, the pairs are not bosons, so that this condensate is not simply 
a Bose-Einstein condensate. 

• The BCS theory predicts an energy required to break a pair. 

• The time-dependent BCS theory also predicts collective modes for the motion of 
the pairs, like sound waves, associated to the famous Anderson-Bogoliubov phonon 
branch [M] [SS] . In a trapped system, the time dependent BCS theory predicts the 
equivalent of superfluid hydrodynamic modes. 

In what follows, we shall use the lattice model of subscction l2'6l In second quantized 
form we therefore take the grand-canonical Hamiltonian 
(169) 

The external potential U(j) may be accompanied by a rotational term if one wishes 
to study vortices. Here the a^r^^'s obey the usual discrete anticommutation relations 
Eqs. (|6l7p . The field operators ipai'^) obey anticommutation relations mimicking the 
ones in continuous space in the limit I — > 0: 

(170) {Mr)Al{r')} = r^5pp,5„,,. 

As discussed in subsection 12 61 this lattice model is very close to the usual Hubbard 
model of condensed matter physics. What is unusual here is that we take a quadratic 
dispersion relation, and more important the model makes sense (to describe a real atomic 
gas with continuous positions) in the low filling factor limit, an unusual limit in solid 
state physics. We note that the original Hubbard model (with the usual cosine dispersion 
relation Eq. (jl6ip and with a filling factor of the order of unity) may be realized in a real 
experiment with atoms in an optical lattice [56l [57l [58l [59] . 



Basic theory tools for degenerate Fermi gases 



37 



3'1. The BCS ansatz. — 

3'1.1. A coherent state of pairs. Let us recall the Glauber coherent state of a bosonic 
field, 

(171) |a) =e-l"l'/2e"''+|0), 

where a is a complex number and creates a boson in some normalized one-body 
wavefunction <p{f), 



(172) a^ = J2l'Hr)^ 



By analogy, the BCS theory [60] , which is a variational theory, takes as a trial state 
a coherent state of pairs: 

(173) |*Bcs) =AAe^^>) 

where Af is a normalization factor, 7 is a complex number and creates two particles 
in the normalized pair wave function (/)(ri, ^2): 

(174) = E I'^i^i^^^) 4'\{ri)i'\{r2). 



Note that in general C and do not obey bosonic commutation relations. For simplicity, 
we omit terms creating two particles in the same spin state. This shall be sufficient for 
our purposes as we restrict here to the case of balanced spin state populations. 

3'1.2. A more convenient form from the Schmidt decomposition. To briefly review the 
main properties of this ansatz, it is convenient to introduce the Schmidt decomposition 
of which is, in the physics of entanglement, routinely applied to the state vector of 
two arbitrary quantum particles: 

(175) 710)-^^"!^")®!^") 



where the coefficients are real numbers, the set of \Aa) is an orthonormal basis here 
for a spinless particle, and the set of \Ba) is also an orthonormal basis for a spinless 
particle. Note that in general these two basis are distinct. As |0) is normalized to unity, 
one has = |7p. Then the pair creation operator takes the simple expression 
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where c^^^ is the creation operator of a particle in the state \ Aa)\ t) and cj^j^ is the creation 
operator of a particle in the state \Ba)\ i). Note that these creation operators and the 
associated annihilation operators obey the usual fcrmionic anticommutation relations. 

The key advantage of the Schmidt decomposition is to allow a rewriting of the BCS 
state, more familiar to the reader: 



(177) 



ll-Bcs) = A/" 



n 



1 + racjj-|-c|j^ 



|0). 



To obtain this expression, we have used the fact that any binary product cjj,| c^^ commutes 
with any c^^Cp^, and the fact that the series expansion of ex-p{'jcl^^cl^^) terminates to 
first order in 7 since (c^)^ = for fcrmions. The form Eq.((T77|) shows that one can 
consider each pair of modes {a t, « 1} as independent, since each factor in the product 
commutes with any other factor. The normalization factor is then readily calculated, 
N = Ha l/\/l + Tq- It can be absorbed in the following rewriting, that may be the one 
directly familiar to the reader: 



(178) 



I* 



BCS/ — 



|0), 



where Ua and Va are the real numbers 

1 



(179) 



Ua = 



and 



Note that they satisfy U^ + V^ = 1, and the minus sign in Va was introduced to ensure 
consistency with coming notations. 

3 '1.3. As a squeezed vacuum: Wick's theorem applies. A third interesting rewriting 
of the BCS state can be obtained from the identity 



(180) 



s^(^^"^-"'')|0) = cos 0|O) + sine 6tct|o) 



where is a real number, b and c are two fcrmionic annihilation operators with standard 
anticommutation relations I f^l Then I^'bcs) — '5'|0) where the unitary operator is 



(181) 



(^) This identity can be proved by a direct expansion of the exponential in powers of 0, or by 
obtaining a differential equation satisfied by the left hand side considered as a function of 6. 
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and the angles 9a are such that 

(182) Ua = cofiOa Va = -sin6'a. 

For bosons, S would be interpreted as a squeezing operator [61]. The BCS state is 
therefore the equivalent for fermions of the squeezed vacuum for bosons. 

Calculating the expectation value in the BCS state of a product of operators c and 
is therefore equivalent to calculating the expectation value in the vacuum state of 
the product of the transformed operators S'^cS and S'^& S. These transformed operators 
have a linear expression in terms of the original c and c^, since one has 

(183) 5'''' Cats' = UaCa^ - ^af^ai 

(184) S'^(^aiS = yo.Cax^Vacl^. 

The linearity of these transformations is evident since S can be formally considered as 
the time evolution operator for a quadratic Hamiltonian, which generates indeed linear 
equations of motion of the creation and annihilation operators. As a consequence, Wick's 
theorem can be applied to calculate expectation values in the BCS state, since it applies 
for the vacuum. 

3'1.4. Some basic properties of the BCS state. Before moving to the energy mini- 
mization within the BCS ansatz, we calculate some interesting quantities. Since Wick's 
theorem applies, the expectation value of any quantity is a function of the only non-zero 
two-point averages: 

(185) (cl^car) = (cl^c^i) = 

^ ~^ ^ OL 

(186) (CciCat) = '\^L{^a\) = -I , p2 ' 

where the expectation values are taken in the BCS state. 

For the total number of particles in the BCS state, we obtain the mean value and the 
variance 



op2 

(187) (N) = E YTT^ 

(188) VariV^^ 

One then immediately obtains the inequality 

(189) Var7V<2(iV), 



40 



YVAN Castin 



an inequality that becomes an equality in the limit where all Fq, <C 1. This shows that, in 
the large N limit, the fact that the BCS state has not a well defined number of particles is 
not a problem for most practical purposes, since the relative fluctuations are 0{1/^/N). 
For the expectation value of the commutator of C and in the BCS state wc obtain 

(190) ([C, C^]) = 1 - ^° , 

2-ia a 

We see that, in the limit where all F^ <C 1, this expectation value is close to unity. In this 
limit, one may consider that creates a bosonic pair: the BCS state becomes a Glauber 
coherent state for this bosonic pair, one obtains Poissonian fluctuations in the number 
N/2 of pairs, which explains the upper limit of Eq. ()189p . This also shows that the BCS 
ansatz will have no difficulty to predict the formation of a Bose-Einstein condensate of 
dimcrs in the 0^0+ limit [52] ■ 

A last interesting property is that the vacuum of an arbitrary quadratic Hamiltonian 
(quadratic in the quantum field) is actually a BCS state. We defer the proof of this 
statement to ^3 31 A similar property is that the vacuum of a set of arbitrary operators 
ba obeying anticommutation relations is a BCS state: in other words, the BCS state is 
a quasi-particlc vacuum. We refer to [63| for a detailed description of this aspect. 

3'2. Energy minimization within the BCS family. - We now proceed with the mini- 
mization of the energy of the lattice model Hamiltonian within the family of not neces- 
sarily normalized BCS states. We define 

(191) ^[^j^(*BCs|i?|*BCS> 



(*BCS|*BCS) 



where the BCS state I^'bcs) is of the form Eq. (|173[) parametrized by the unnormalized 
pair wavcfunction $ = 70 and by the factor f\f now considered as an independent variable. 
If one performs a variation of $ around a minimizer of £'[$], $ = $o + 6^, this induces 
a variation i5|5'bcs) of the BCS state around I^'bcs)' 



(192) 5|*BCS> = J2 l'SHri,r2)i'\{ri)i^lir2)\Kcs) = SX\^b 



BCS/ 



The corresponding first order variation of the energy function E[^] has to vanish for all 
possible values of S^, 

(193) 5£;=((5(*Bcs|)(i^-^[*o])|*Bcs)+c-c. = V 5$, 

where we assumed that the minimizer I^'bcs) normalized to unity. 

The second step is to introduce the quadratic Hamiltonian deduced from the full 
Hamiltonian H by performing incomplete Wick's contractions. This recipe has to be ap- 
plied to the interaction term only, since the other terms of the Hamiltonian are quadratic. 
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The quartic on-site interaction is replaced by the quadratic expression: 

(194) +^|^Tffo(^I^i)+ !<->!, 

where the expectation values are taken in the state l^'gcs)- '^^'^ ^^^^ ^i-"^*^ ^^^'^ righthand 
side consists in a pairing term, that creates/annihilates a pair of particles; it involves the 
following pairing field, also called 'gap' for historical reasons that will become clear: 

(195) Air) = go{i^iif)^^ir)). 

The second line of the righthand side is simply the Hartree mean field term: each spin-up 
particle sees a mean field potential which is go times the density of spin-down particle. 
Note that, in the lattice model, this mean field term does not diverge at the unitary 
limit, since qq remains bounded in this limit, see Eq. (jl55p . 

The quadratic Hamiltonian associated by Wick's contractions to the full Hamiltonian 
is the so-called BCS Hamiltonian. and it has the structure 



(196) H = ^ /3 i^Ur)hr-r-'Mr') + ^ [^j(r )ViI(r ) A(f ) + h.^ 



iadj- 



Here the matrix h is the sum of the one-body part of the Hamiltonian H and of the 
Hartree mean field terms. 



(197) hp^p, = [kin]- + U{r) + goiV-j (^)V^t(^)) 



where [kin] is the matrix representing the kinetic energy operator. We have assumed 
for simplicity that the one-body Hamiltonian is spin independent and that the two spin 
density profiles are identical, so that the matrix h does not depend on the spin state. 
Finally, the additive constant E'adj is adjusted to the value 

(198) i^adj = -goY.^' [(^t^t)' + K^xV-t)!' 

r 

to ensure that the mean energy of Ti and of H coincide in the minimizer I^bcs)' 

(199) {H) = (H). 

Using Wick's theorem and the fact that the operator 6X is a quadratic function of 
the quantum field, one can then show the marvelous property 

(200) iS(^BCs\)H\Kcs) = ('^(*Bcs|)H|vI/Ocs), 



whatever the variation (5|^'bcs) of the BCS state within the BCS family. 
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There is the foUowing key consequence. Considering the quadratic Hamihonian Ti. 
as a given Hamiltonian, it turns out that the BCS state I'I'bcs)' obeying Eq. (|193|l . also 
obeys the condition to be a stationary state of the energy functional 

(201) Em = (^BCS|H|V&BCS) 

(^'bcsI^'bcs) 

noting that one has i?[(f>o] = f [$o]- Condition Ea. p93p is thus equivalent to the condition 

(202) <5£ = (<5(*Bcs|)(W-f[$o])|*Bcs>+c-c. -0 V 5$, 

That is, to minimize within the BCS family, we shall take the BCS state that 

minimizes the energy of Ti. 

We are left with the study of a quadratic Hamiltonian. As reminded to the reader in 
ij3'31 the quadratic Hamiltonian Ti may be diagonalized by a canonical transformation, 
the so-called Bogoliubov transformation. This allows to show that the ground state 
I ground) of TL is in fact a BCS state (see ^3'3p . and to calculate the expectation values 
(i/'J.'0CT)ground and (V'xV't) ground in |ground). Since 1^*1303) ^^"^ |ground) actually coincide, 
one is left with the self-consistency conditions: 

(203) (^i^/ia>grou„d = (V^i^.) 

(204) (^i^T)grou„d = (V^X'/'t)- 

Since the expectation values in the ground state of Ti. are non-linear functionals of the 
coefficients {'ip'^'ijjc) and ('0xV't) of this constitutes a non-linear self-consistent problem, 
the so-called BCS equations for the density and the gap function. 

Another crucial consequence is to realize that one may write the self-consistency 
conditions taking an excited eigenstate of the quadratic BCS Hamiltonian Ti rather than 
the ground state in this case, the variational BCS calculation is used to predict 

excited states of the gas! This shows that BCS theory immediately predicts elementary 
excitations of the gas. As we shall see in the thermodynamic limit, these excitations 
correspond to pair breaking and are characterized by an energy spectrum with a gap. 

3'3. Reminder on diagonalization of quadratic Hamiltonians . - We consider here a 
quadratic Hamiltonian of the form Ea. p96p . where, at this stage, h is an arbitrary (spin 
independent) hermitian matrix and A(r ) an arbitrary complex field. We wish to put 
this Hamiltonian in canonical form by performing a Bogoliubov transformation. 

The key property of a quadratic Hamiltonian is that it leads in Heisenberg picture to 
linear equations of motion. Collecting the field operators and their hermitian conjugate 
at all points of the lattice in big vectors, we get 



(205) ihdt 
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where —a stands for the spin component opposite to cr and where we have introduced 
two matrices 

h diag(A) \ h diag(-A) \ 

diag(A*) -h* ) \ diag(-A*) -h* ) ' 

These two matrices are hermitian: they can be diagonaUzed in an orthonormal ba- 
sis. Furthermore they obey the following symmetry property: if the vector (w, v) is an 
eigenstate of with the energy e, then 

• the vector (— w*, u*) is an eigenvector of ^ with the eigenvalue — e, 

• the vector (u, ~v) is an eigenstate of L^, with the eigenvalue e, 

• the vector [v* , u*) is an eigenvector of , with the eigenvalue — e. 

Assuming for simplicity that all the cigencncrgies are non-zero, we can collect the 
eigenvectors in pairs of opposite eigenenergies. Expanding in the eigenbasis 

of L'^ , and expressing the negative energy modes of and all the eigenmodes of in 
terms of the positive energy modes of L^, noted as {ua,Va) , we obtain 

(207) i>^{f) 

(208) Vi^(f) 

where the positive energy eigenvectors of arc normalized in a way mimicking the 
continuous space limit, with Dirac's notation {f\ua) = Ua{r): 

(209) (ua|u/3) + {va\vp) = ^ < (r )u/3 (t^) < (r )i;^ (r ) = S^p, 

f 

The coefficients in the expansion of the field operators are themselves operators that are 
easy to express using the orthonormal nature of the eigenbasis of L"^: 

(210) = Z'E<(^^)^t(0 + <(^^K'I(^') 

r 

(211) 6«i =/3^<(f)^x(r-)-<(r)^|(f). 

f 

Using Ea. (|209p again, one can check that the 6's and their hermitian conjugates obey 
the standard fermionic anticommutation relations. 

The last step is to express the Hamiltonian in normal form, using the 6q,'s. One first 
realizes that the quadratic Hamiltonian can be expressed as a quadratic form defined by 



(206) = 



= ^foaTWa(^) - &l|<(r) 

a 
a 
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the matrices L"^ , 

(212) 7i = Tr/z + i?,dj + ^/^X^(V^-V^-)^"( /t" ) 

where the extra constant term, given by the trace of the matrix h, originates from the 
anticommutator of ip^ with ip^. Inserting the modal decomposition Eqs. (|207l208p and 
using the fact that (mq,Wq) is an eigenvector of L^, etc, we obtain the canonical form: 

(213) 71 = Eo + Y,e^bi,ba,a. 

a,cT 

where we recall that all the Cq, > by construction. The immediate expression for Eq is 
Eq = Tr /i + i?adj — ^"i involves two quantities that are not well behaved in the zero 
lattice spacing limit. Expressing the matrix as a sum of ieo, times dyadics involving 
the eigenvectors as ket and bra, and projecting over the upper left block leads to 

(214) hf^^p, = E K(r )<(r') - <(r )i;„(r ')] , 

a 

where we used the orthonormal nature of the cigcnbasis of , as defined in Eq. (j209p . 
Taking the trace of this expression in the spatial basis of the lattice leads to the more 
convenient form for Eq: 

(215) E^=E,i;-2Y, 

a 

where E^adj is given by Eq. p98p . The canonical form Eq. (|213p . with all the Cq, > by 
construction, immediately shows that the ground state |^'o) of Ti. is the vacuum of all 
the 6q.'s, with an eigenenergy Eq. Please remember that we are dealing here with the 
grand canonical Hamiltonian; this is why £^0 < for the ideal gas go = 0. 

To show that this ground state is a BCS state, we take for |^'o) an ansatz of the BCS 
type Eq. (|173p to solve the equations 

(216) &a<x|*o)-0 

for all mode index a and spin component a =|, j. The commutator of a baa with the 
pair creation operator is a linear combination of the ?/;t's so it commutes with C^. As 
a consequence, 

(217) [6„,,e^^'] =7e^^'[&a.,C+]. 
Since exp(7C^) is an invertible operator, Eq. (|216p reduces to 

(218) {b^a+l[bc.a,C^)\0) =0. 
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From the expressions Eas. (|210l211[) we calculate the commutator with and obtain the 
linear system 

(219) <(r) +75]/W,r)<(r') = 

r' 

(220) (r ) + 7 E ^'X (^^') = 

r' 

The solution may be written in matrix form as follows. Since the lattice has discrete 
positions, wc consider 0(r, r ') as the element of a square matrix [0] on row r and column 
r '. Similarly we form the square matrix \U] (respectively \y]) such that the element on 
row r and column a is Ua{r) (respectively Va{f)). Then one has the explicit writing 

(221) ^[cf\=^'[<p]^-l-^[U\^)-\V]\ 

where *M is the transposed matrix of a matrix M . Since the set of all the eigenvectors 
of forms an orthonormal basis, one has I'^iU'^U + V'^V) = Id and *UV = *VU ^ where 
Id is the identity matrix. This ensures that the solution [0] is indeed a symmetric matrix 

(222) (/'(ri,r2) = 0(7^2, ri). 

This symmetry property could be expected from the fact that the quadratic Hamiltonian 
Ti is invariant by an arbitrary rotation R of the spin degree of freedom, and that its 
ground state is not degenerate since we have assumed Cq strictly positive. The condition 
that = RC^ for all spin rotations indeed imposes that </> is a symmetric function 
of ri and r2 ■ 

Using Wick's theorem one can easily calculate expectation values of observablcs in 
I'I'o) using the modal decompositions Eqs. (|207l208p . Two important examples are the 
mean density and the anomalous average 

(223) p^if)=p^if)=Y,Kir)\\ 

a 

(224) {Mr)i'lir)) = - J] «a(r )<(f ). 

a 

3'4. Summary of EC'S results for the homogeneous system. - 

3'4.1. Gap equation in the thermodynamical limit. In the case of a spatially homo- 
geneous system (cubic box with periodic boundary conditions) explicit predictions are 
easily extracted from the BCS theory. For a spatially uniform gap parameter A, assumed 
to be real non negative, and for a uniform density profile p| = pj^ the spectral decompo- 
sition of is readily obtained: from the translational invariance one seeks eigenvectors 
of the form 



(225) 



(«j(f),«,,(f)) = L-3/V'^-(C/,,y,). 
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Restricting to the positive part of the spectrum, one finds the eigenenergies 

1/2 



(226) 



2 1,2 



2m 



2z,2 



h'k 



2m 



fi + iA 



where jl ~ /j — ,9oPt chemical potential minus the Hartree mean field term. The 

amplitudes Uk,Vk, chosen to be real, are normalized as C/| + = 1 and are given by 



(227) 



To make the link with the quantities Fq of section [3'1.2i we take for \Aa) = |v4g) the 
plane wave of wavevector k and for \Ba) = |-B^r) the plane wave of wavevector —k. This 
gives 



(228) 

and a pair wavefunction [65] 



— - A* + efc 



(229) 



fee 



ik-{ri —r2) 



keV 



Remember that the regime where all the <C 1 corresponds to the case where the BCS 
state is a condensate of bosonic pairs. 

Taking the real and imaginary parts of Eq. (|227p . and going directly to the thermo- 
dynamic limit, one obtains from Eqs. (|223|224p 



(230) 
(231) 



d?k 
(2^ 



1 - 



A 



-.90 



(Pk A 
P (27r)3 2^' 



where p is the total density. This makes explicit the self-consistent character of the 
conditions Eqs. (|203l204p . The last equation Eq. (|23ip is called the gap equation. 

3'4.2. In the limit of a vanishing lattice spacing. An important point is to discuss the 
dependence of these predictions with the lattice spacing I, since the relevant regime of 
a continuous gas is described by the lattice model only in the I — > limit. The integral 
giving p converges at large fc, since the integrand is 0(l/fc''), so that one may directly set 
Z in Ea. (|230p . On the contrary the integral in Ea. (|23ip has an ultraviolet divergence 
when integrated over the whole momentum space. However tends to zero in the / 
limit, compensating exactly this divergence. Dividing Ea. (|23ip by and expressing 1/ go 
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from Eq. (jl54p leads to a more convenient expression for the gap equation, in which one 
may take directly I 0: 



(232) — = A 



9 Jv (27I-) 



m 1 



Furthermore, we note that the BCS theory becomes simpler when / 0: since the 
coupling constant go vanishes in this limit, the Hartree mean field term also vanishes, so 
that fl ^ fi and the matrix h in the quadratic Hamiltonian Ti. reduces to the matrix of 
the kinetic energy minus fi. The a priori unknown coefhcients of Ti. thus reduce to the 
anomalous average (ip^tp-^), and Eq. (|204p remains as the only self-consistent equation, the 
other one Eq. (|203[) giving explicitly the density as a function of /i and of the anomalous 
average. 

In general, one has in practice to use the BCS theory in this I — s- limit, to be in the 
continuous gas regime. As we shall see in section [3"4.4[ a notable exception is the weakly 
attractive regime kpa — > 0~. 

3'4.3. BCS prediction for an energy gap. As sketched in the last paragraph of 
section 13 2[ the excited states of the BCS Hamiltonian H correspond to a variational 
BCS approximation for excited states of the full Hamiltonian of the gas. Here we come 
back to this statement, and show how to derive one of the most celebrated predictions 
of BCS theory, the minimum energy required to break a pair. 

First we revisit briefly section 13 21 in the process of energy minimization within the 
BCS family of states, we ended up looking for a BCS state that is the ground state of a 
quadratic Hamiltonian of the form Ea. (|196p . Taking here for simplicity the continuous 
limit ^ ^ 0, for a spatially homogeneous solution, we can omit the Hartree mean field 
term and perform a variational calculation that immediately restricts to BCS states being 
the ground state of the following Hamiltonian, 

(233) = 5] [Ek - ^i) a\^a^^^ + l^Y. H (^^) + ^'C' 



where = h^k^/2m and the real non-negative number A is now the only variational 
parameter. The expectation value of the full Hamiltonian H in the ground state Ivt*^) of 
Ti.\ is readily evaluated from section [3^ and from Eas. (|226l227p . 

(234) iH)i ,) (i - ""r.'ai ) + -^"^^ i^^^^^^^'^T' 

with the anomalous average (V'j.V't)^ ~ —L^"^ J2kev ^/ i'^\Ek ^ M + In these expres- 
sions, the appex g stands for "ground state" . It remains to require that the first order 
derivative of with respect to A vanishes, which reduces, after explicit calculation of 
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the derivative, to the foUowing equation satisfied by the optimal value Xg of A: 

(235) 5o(^iV'T)t = A,. 

One has recovered, in a slightly different way, the self-consistency condition Eq. (j204|l . 
and one finds that Xg satisfies the gap equation (with jl = ji). This optimal Xg was called 
earlier A. 

We can now perform a similar variational calculation, taking as a variational ansatz 
an excited state of T-Lx . To preserve the symmetry between the two spin components t 
and [, we consider the excited state 

(236) \^l)=b\^bl.^\^i), 

where 9 is a fixed wavevector. Two Bogoliubov modes of the Hamiltonian Tix are now 
in the excited state, with occupation number one, rather than being in the vacuum state 
of the corresponding 6's. The mean energy of H in this excited state is simply, in the 
thermodynamic limit, 

(237) ^ m + 2(i^,- M)- + 2Ago(j/,^T)-. ^ ^(^-3). 

\Eg-fl + tX\ 

It differs from (i?)^ by a term 0(1), whereas the full energy is 0{L^), so that the optimal 
value Ae of A for the minimization of differs from Ag = A by a small term 5X. We 
can then simply expand Eq. (j237p in powers of SX. The key point now is that the ground 
state energy varies only to second order in dX since it is minimal in A; since the second 
order derivative of (i^)^ is O(i^), 6X is 0{1/L'^) and the contribution to the total energy 
of the corresponding second order variation of (i/)^ is 0{1/L^) and negligible. In the 
remaining terms of Ea. (|237p we can simply set A = A. Using Eq. (|235p . we are left with 

(238) {H)l^{Hy^ + 2e, + Oil/L^) 

where Cq is the BCS spectrum Ea. (|226p . 

To summarize, one considers excited states of Ti. as variational states, with a small 
number of Bogoliubov excitations. One should then in principle solve again the self- 
consistency conditions Eqs. (|203l204p . but in the thermodynamic limit the excitation of a 
few Bogoliubov modes has a small effect on the density and gap parameter, an effect that 
we have shown to be negligible on the energy. One can then directly consider that the 
elementary excitations of the quadratic Hamiltonian Ti. associated to the ground energy 
BCS state actually give the energy of the elementary excitations of BCS theory [M] . 

When A > 0, we see that the minimal value of the BCS spectrum with respect 
to k is non-zero: it has a gap -Egap ~ niin^ej;. When /i > 0, which contains the usual 
regime of the BCS theory, the regime of condensation of Cooper pairs (see below), the 
gap is iJgap = A, hence the name of A. When /2 < 0, which contains the regime of a 
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Bose-Einstein condensate of dimers, the gap is Sgap = (A^ + A^) . Finally, we get 
from BCS theory a prediction for the miminal energy 2iJgap required to break a pair, 
that is to get from a condensate of N/2 pairs, a condensate of N/2 — 1 pairs and two 
unpaired atoms |67| . 

3'4.4. BCS predictions in limiting cases. We shall not discuss the full solution of 
the BCS equations but we briefly review simple limiting cases. We introduce the Fermi 
momentum kp of the ideal Fermi gas, such that k'p/i/oi:'^) = = p^. 
Limit kpa 0^: the gap parameter A tends exponentially to zero in this limit, as we 
shall see. For such a small value of A, one may set A = in Eq. (j230p . to obtain that p. 
is the Fermi energy of the ideal gas h^kp/2m, so that 

(239) (67rV)'/'+ffoPT- 

2m 

If one takes the mathematical limit I — *■ the Hartree mean field term disappears since 
go ^ in this limit, as already mentioned. However in the present weakly attractive 
limit one can choose |a| ^ / ^ fc^^ so that the continuous space limit and the zero-range 
condition Ea. (|116p for the lattice model are obtained while go — g, see Ea. (|155p . The 
on-site interaction potential is then treatable in the Born approximation, which makes 
the BCS approach more accurate. One then indeed recovers the first term gp-^ of a 
systematic expansion of p in powers of kpa [29] . Now turning to the gap equation, one 
finds a gap parameter [45| 

(240) A~ 8e-2/ie-'^/(2''^^l''l^ 

Since p = h^kp/2m > in this weakly attractive limit, \Tk\ can assume extremely large 
values, for h'^k^/2m < p: the pairs that are condensed are not bosons. Note the non- 
analytic dependence of the gap on the small parameter fcF|a|j which indicates that the 
BCS state in the thermodynamic limit cannot be obtained by a perturbative treatment 
of the interaction potential. This non-analytic dependence can be readily seen from 
Eq. (|232p . whose integrand diverges in A; = ^/2mp,/h, in the limit A — ^ 0. Replacing the 
k integration variable by the ideal gas mode energy E = h^k'^/2m, and approximating 
the ideal gas density of modes by a constant in the energy interval of half-width (5 <C /i 
around p, one gets a contribution 

(241) r , "^"-"'f 

The same technique applied over e.g. the interval E E [0, 4p] leads to Eq. (|240p . when one 
also takes into account (to lowest order in A) the difference between the exact integral 
in the gap equation and the approximate one. 

Limit kpa 0+: the coefficients F^ tend uniformly to zero in this limit, because 
p < and A <C \p\, so that the pair creation operator obeys approximately bosonic 
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commutation relations. This was expected physically, since the ground state of the gas is a 
condensate of almost bosonic dimers. Since two atoms in a dimer are at a relative distance 
a, one should ensure ^ ^ a in the lattice model, so that one takes the mathematical 
limit I — > 0, which implies /i — > ^. Let us check that the BCS theory correctly predicts 
a condensate of such dimers. We first simplify the gap equation Eq. (|232p by using 
ek — — /i + h^k^/2m. After division by A, we get an equation for the chemical potential 



(242) 

which leads to 



1 r (fk 



9 J (27r) 



1 



(243) ■ 



- 2u ?i^A;2 

m 1^ 



This is minus half the binding energy of a dimer, exactly what was expected (keep in 
mind that ^iN = A^moi^moi where iVmoi is the total number of dimers and is equal to 
N /2, so that the molecular chemical potential /imoi is twice the atomic one). The next 
step is to expand the integrand of Ea. (|230p to leading order in A to calculate the gap 
parameter [45] : 

(244) A^J^i^k,afl^^«2L 

Note that in this molecular BEC regime, A is not proportional to the energy required to 
break a pair. One sees from Eq. (|226p that the gap in e/c is ~ since /i is negative and 
much larger in absolute value than A. The energy to break a pair is then 2|^|, which is 
indeed the binding energy of a dimer. Finally, by performing the Fourier transform in 
EQ. (|229p . approximating Vk/Uk to leading order in A, one obtains the pair wavefunction 

(245) ^{rl - ~ j^M\n - r1|) 

where 0o is the wavefunction of the bound state of two atoms given by Eq. (|146| ). and 
where we used 7^ = J2ki^k/Uk)'^ . 

Limit fecial = +00: the numerical solution of the gap equation Eq. (j232p and of the 
density equation Eq. (|230[) in the limit I gives the BCS estimate of the numerical 
coefficient 77 of Ea. (|168p for the groimd branch. This estimate is an upper bound [55] : 

(246) Tj < mcs = 0.5906 . . . 

A much better estimate was obtained by approximate (fixed node) Monte Carlo cal- 
culations [351 136] . 77 ~ 0.4, a value confirmed by recent exact quantum Monte Carlo 
calculations [33] . Early measurements of 77 were in contradiction with these theoretical 
values [37] , but more precise measurement performed in Innsbruck [38] , in Paris [39] and 
at Rice University [43] are consistent with them. 
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3'5. Derivation of superfluid hydrodynamic equations from BCS theory. ~ A key point 
of the BCS theory is to predict a spectrum of elementary excitations having a gap, see 
Ea. (j226[ ). One could then be tempted to infer strong predictions on the thermodynamic 
properties, e.g. that the entropy obeys an activation law ©(e^^sap/'^B'^") at low tempera- 
ture, where i?gap is the energy gap. However, in the BEC limit, where the gas is simply 
a Bose-Einstein condensate of bosonic dimers at low temperature, this prediction of an 
activation law is obviously wrong, since one knows that the relevant excitation spectrum 
is the Bogoliubov spectrum, which has no gap. 

In reality, the spatially homogeneous gas, whatever the considered regime (BEC or 
BCS), is expected to have a branch of collective excitations with no gap, behaving as 
sound waves at low momenta. These collective excitations correspond to coherent os- 
cillations of the density of the pairs, which are not gapped, to be distinguished from 
the pair breaking excitations which are gapped. The key point briefly addressed in this 
section is that the time- dependent BCS theory contains such a branch of collective, non- 
pair-breaking oscillations |69| . This is similar to what happens for the weakly interacting 
Bose gas: the Bogoliubov spectrum is obtained from a linearization of the time dependent 
Gross-Pitaevskii equation. 

Here, rather than performing an exact linearization of the time-dependent BCS equa- 
tions, which leads to the so-called RPA approach [531 EOl ED EH ESI El] j we go through 
a sequence of simple approximations allowing one to derive superfluid hydrodynamic 
equations from the time-dependent BCS theory. This has several advantages: it is more 
physical, it easily applies to harmonic trapping [75], and it is applicable to the non- 
linear time dependent regime, to study analytically low-energy collective excitations of a 
trapped superfluid Fermi gas |7S1 [77] , but also the expansion of the gas after the trap was 
switched off [75] and the response of the gas to a rotation of the harmonic trap [TS] [50] 
even in the non- linear regime [80] . 

Coming back to thermodynamical aspects, one may fear that the straightforward finite 
temperature extension of the BCS theory, with a variational density operator which is 
Gaussian in the field operators |63| . is not able to calculate the critical temperature Tc 
with a good accuracy, because it does not correctly include the collective excitations. It 
turns out that this simple finite temperature BCS theory correctly gives Tc in the BCS 
limit kpa 0^ within a numerical factor [81j . but is indeed totally wrong in the BEC 
limit kpa 0^ where it does not reproduce at all Einstein's prediction for the critical 
temperature of the ideal Bose gas. Refinements of the BCS theory have been developed 
to recover Einstein's prediction and to obtain a calculation of Tc approximately valid 
over the whole range of values of kpa [HI H5] . 

3'5.1. Time dependent BCS theory. This is a direct generalization of the static case 
of section I3'2I The exact A^-body Schrodinger equation can be obtained from extremal- 
isation over the time-dependent state vector \^{t)) of the following action. 



(247) 



S 



dt<- 



(*WI^I*W)-c.c. 



{^{t)\H{t)\^{t)) 
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for fixed initial and final values of the state vector. For an arbitrary 

variation \S'i>{t)) of the ket \'i>{t)), subject to the condition 



(248) 



iJM'(to) = mm = 0, 



we calculate the first order variation of the action, 



(249) 



SS 



dt \ (,5^'(t)| 



c.c. 



We have integrated by parts over t to get rid of {d/dt)\S'i') and we have used Ea. (|248p 
to show that the fully integrated term vanishes. The condition that SS vanishes for all 
ket variations indeed leads to Schrodinger's equation. 

For time dependent BCS theory, one forces the ket to be of the BCS form Eq. (|173p , 
I'J') = I^BCs)- The variation of the ket is performed within the BCS family, by a 
variation of $ = 70 and of A/", so that \S'9) ~ <^(|^'bcs))- Wc can then introduce the now 
time dependent quadratic Hamiltonian Ti.{t) constructed in section [3 ' 21 such that, for all 
variations within BCS family. 



(250) (<S(*Bcs(i)l)^WI*Bcs(i)) = ('5(*Bcs(t)l)H(t)|*Bcs(t)>. 

One then sees that SS in Eq. (|249p identically vanishes if the ket evolves with the quadratic 
Hamiltonian. 

(251) »?i^|*Bcs(i)> =WWI*Bcs(i)>. 

It remains to check that the ket evolving this way remains of the BCS form. To this 
end, one shows, with the reasoning having led to Eq. (|217|) . that ipcr\'^BCs{t)) is equal to 
a linear combination of the ipl,{r) acting on |^'BCs(i))- One then sees that H|^bcs) is 
of the form "constant plus polynomial of degree exactly two in ip''" acting on |vE'bcs(0)j 
which can be reproduced by an appropriate time dependence of $. 

In practice, the equation of motion for $ is not useful. One rather moves to the 
Heisenberg picture with respect to the time dependent Hamiltonian Ti.{t), as we did in 
section I3'3I Since this Hamiltonian is quadratic, the equations of motion for the fields 
are linear, of the form Ea. (|205p . where is now time dependent. Since these equations 
are linear, we can solve them by evolving in Eas. (|207l208[) the mode functions {ua,Va) 
while keeping constant the quantum cocfhcients baa where a =t or J,: 

(252) ma,(-).Ltw(-). 

These equations on the mode functions are effectively non-linear since coefficients of , 
involving the mean density and the anomalous average {tpi'ip'^), depend on the 
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mode functions; since the baa are constants of motion, this dependence is stiU given by 
Eqs. (|223|224|l . now taking the time dependent u^'s and w^'s. 

3'5.2. Semi-classical approximation. The physical situation that we have in mind here 
is a gas in the BCS regime (a < 0, /i > 0) in a time dependent harmonic trap. The trap 
may be rotating, with an angular velocity fl{t), in which case one moves to the rotating 
frame to eliminate the time dependence of the trapping potential due to rotation; this 
introduces an additional term —n{t) ■ L in the one-body Hamiltonian, where L is the 
angular momentum operator of a particle. In this case one may expect that quantum 
vortices form j82| for a high enough rotation frequency; the semi-classical approximation 
that we present is however restricted to a vortex-free gas. 

We wish to treat the equation of motion of (ua,Va) in the semi-classical approxima- 
tion. The general validity condition of a semi-classical approximation is that the applied 
potentials vary slowly over the coherence lengths of the gas, a coherence length being 
the typical width of a correlation function of the gas. 

A first correlation function here is Using the homogeneous gas expres- 

sion of this function, one sees that it is the Fourier transform of the function of width 
the Fermi momentum, Afc ~ kp, where ?i^fc|,/2m = /i. The associated coherence length 
is 1/Afc = kp^. 

A second correlation function is ('^'))- For the homogeneous gas this is the 

Fourier transform of the function UkVk, which according to Eq. (|227p has a momentum 
width Afc = m\A.\/TT?'kp. The associated coherence length 1/Afc corresponds in BCS 
theory to the length of a Cooper pair, 

(253) Zbcs 



toIAI 



Since A < /i in the BCS regime, we keep /bcs as the relevant coherence length. 

A first typical length scale of variation of the matrix elements in Eq. (|252p originates 
from the position dependence of |A|: in the absence of rotation and for an isotropic 
trap, this is expected to be the Thomas-Fermi radius i?TF of the gas, defined &s ji ~ 
where ui is the atomic oscillation frequency. This assumes that the scale 
of variation of the modulus of the gap is the same as the one of the density, a point 
confirmed in the adiabatic approximation to come. The condition that the mean field 
Hartree term and the harmonic potential have a weak relative variation over ^bcs for 
typical values of the position also leads to the condition 

(254) Ibcs < Rtf- 

For an isotropic harmonic trap this is equivalent to the condition 



(255) 



|A| > huj. 
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In the general time dependent case, however, this is not the whole story, since the 
phase of A can also become position dependent. In the rotating case for example, with 
a rotation frequency of the order of to, the phase of A may vary as ~ mujxy/h; this 
introduces a local wavevector mcuRT^p/h kp, making a semi-classical approximation 
hopeless. We eliminate this problem with a gauge transform of the u's and v's |80j : 

(256) S„(r,t) = u^{r,t)e-'^^'^'*^^^ 

(257) ^}„(r,^)EE«„(f,^)e+'^(-~*)/^ 

where S is defined as 

(258) A(f, t) = |A(f, t)\ e2^•5(^^^*)/^ 

The time dependent BCS equations Ea. (|252p are modified, A being replaced by |A| and 
the single particle Hamiltonian h being replaced by the gauge transformed Hamiltonian 

(259) h = e~'^^''he+'^^^ + dtS. 

An efficient frame to perform a semi-classical approximation in a systematic way for 
the evolution of a wave in a potential is to use the Wigner representation [6] of the density 
operator of the wave. Here the wave is represented by the two-component spinor (mq, Vq,) 
so that we introduce the corresponding density operator 



(260) 



0-11 0-12 \ = ^ ( \^a){Ua\ \Ua){Vo 
<5'21 0-22 / ~ ^ V \Va){u 



Note that the matrix elements of a in position space are related to the previously men- 
tioned correlation functions of the gas, up to the gauge transform, which makes the 
discussion consistent: 

(261) \{r\^22r)\ = \{^ir)^^{r'))\ 

(262) \{r\a,,\f')\ = \{Mr)^lir'))\- 

We introduce the Wigner representation of ct [5] , assuming for simplicity a continuous 
position space: 

(263) Wir,p,t) = Wigner{<7} = / -i^(r - x/2\a\r + x/2)e'P•^/^ 

J [2nhy 



Since a is the density operator of a two-component spinor, the Wigner distribution W 
is a two by two matrix. Within this representation, the key quantities of BCS theory. 
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the total density, the modulus of the gap function and the total matter current have the 
following expressions in the general case of a rotating frame: 



(264) 






f d''pW22{r,P,t) 
hv 




(265) 


|A|(f,t) 


= -go / d^pWi2{f,p,t) 
JhV 


(266) 




= P 


v{r,t) - Q{t) X f 





hv 



where f] is the angular velocity of the rotating frame, and the so-called velocity field is 
defined as 



(267) 



v{r, t) 



gradS'(r, t) 



We have introduced here the total matter current J(r, i), that obeys by construction 
(268) dtP + AYV] = Q, 

In the rotating frame, in a many-body state invariant by exchange of the spin states t 
and [, it is very generally given by 



(269) 



(^'T(^Ograd-0t(r,t)) 



(0(f, t) n{t) X f. 



Within BCS theory, these last two relations still hold [53j and lead to Eq. (|266p . 
The semi-classical expansion then consists e.g. in 



(270) 



Wigner{V^(f')a} = [V{r) + '-^dpV • . . .]Wif,p,t), 



where 1^ is a position-dependent potential. The successive terms in this expansion we 
call zeroth order, first order, etc, in the semi-classical approximation. Assuming that the 
momentum derivative of is ~ W/Ap, where Ap is the momentum width of W, we 
recover the fact that the first order term is small as compared to the zeroth order one if 
the variation of V over the coherence length h/Ap, which is of the order of the spatial 
derivative of V times the coherence length, is small as compared to V. 

Here we shall need only the equations of motion of the Wigner distribution W up to 
zeroth order in the semi-classical approximation: 



(271) 



ihdtW{f,p,t)^ Llir,p,t),Wir,p,t) 
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where the two by two matrix Lq is equal to 

(272) Llir^p^t) =( ^^«(^^'*) 1^1^^^'*) ] . 

° ^ |A|(r,i) -iL-+^,ff(f,i) y 

We have introduced the position and time dependent function. 

(273) fi^sir, t) = i-L- U{f, t) - ^mv^{f, t) + mv{f, t) ■ {n{t) x f ) - dtS{f, t), 

that may be called effective chemical potential for reasons that will become clear later. 

At time < = 0, the gas is at zero temperature. By introducing the spectral decompo- 
sition of {t ~ 0) one can then check that 

(274) a{t ^ Q) ^ 0[LHt = 0)] 

where 6{x) is the Hcavisidc function. Since lJ(< = 0) is the classical limit of the operator 
{t — 0), the leading order semi-classical approximation for the corresponding Wigner 
function is, in a standard way, given by 

(275) W{r,p,t^ 0) cj^9[Llif,p,t^ 0)]. 

This implies that the 2x2 matrix W is proportional to a pure spin-1/2 state \ip){tp\ with 



(276) \^ir,p,t^O)) = 



Uo{r,p) 
Vo{r,p) 



where (C/q, Vb) is the eigenvector of L\{r,p^t = 0) of positive energy and normalized to 
unity. At time t, according to the zeroth order evolution Ea. (|27ip . each two by two 
matrix W remains a pure state, with components U and V solving 

(277) ^wA 

V y{f,P,t) J ^ \ V{r,p,t) 

3 '5. 3. Adiabatic approximation. We then introduce the so-called adiabatic approx- 
imation: the vector {U,V), being initially an eigenstate of L\{r,p,t = 0), will be an 
instantaneous eigenvector of Ll{f,p,t) at all later times t. This approximation holds 
under the validity condition of the adiabaticity theorem |83| . discussed for our specific 
case in |80] . generically requiring that the energy difference between the two eigenval- 
ues of lJ (r, p, t) be large enough. As this energy difference can be as small as the gap 
parameter, this imposes a minimal value to the gap, not necessarily coinciding with the 
one of Eq. (j255p . as shown in [80'. In this adiabatic approximation, one can take 

(278) W{r,p,t) ^^-1^9[Ll{f,p,t)] = \ + (f,p, t)) (+(f,p, i)| 
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where | + (r, p, t)), having real components (CAnst, VJnst), is the instantaneous eigenvector 
with positive eigenvalue of the matrix Lq defined in Eq. ()272| l. (?7inst, Vinst) are simply 
the amplitudes on the plane wave exp(ip ■ r/K) of the BCS mode functions of a spatially 
uniform BCS gas of chemical potential ^cff and of gap parameter |A(f, t)|. Using Eq. (|264p 
and Eq. p65p with the approximate Wigner distribution Eq. (|278p . one further finds that 
this fictitious spatially uniform BCS gas is at equilibrium at zero temperature so that 
the zero temperature equations of state may be used, relating the chemical potential to 
the density, /ip [p] , and the gap to the density, Aq [p] , where the functions //q and Aq may 
be calculated from the spatially homogeneous BCS theory. In particular, the equation 
of state relating the chemical potential to the density leads to 



which leads, together with Ea. (|273p . to one of the time dependent hydrodynamic equa- 
tions, the Eulcr-typc one. Also, t/inst and VJnst are even functions of p, so that the integral 
in the right hand side of Eq. (|266p vanishes and Eq. (|268p reduces to the hydrodynamic 
continuity equation in a rotating frame. 



Under the adiabatic approximation, the superfluid hydrodynamic equations are thus 
derived from BCS theory, remarkably without having to postulate the expression of the 
superfluid current in terms of the gradient of the phase of the order parameter. We refer 
to [80| for a discussion of the validity condition of the adiabatic approximation. 

The presence of superfluid hydrodynamic equations allows to conclude in the existence 
of collective modes in the BCS theory, in a transparent way, by a simple linearization 
around steady state. Although we have considered here a trapped system, we note that 
the formalism may also be developed in the spatially homogeneous case. To create a 
sound wave of wavevector q, one may apply an external potential varying with spatial 
harmonics e^"'''', e.g. with a Bragg pulse [HH [HS]. The semi-classical approximation 
holds if this external potential varies slowly at the scale of the Cooper pair size ZbcSj 
that is (//bcs ^ 1- Linearization of hydrodynamic equations in the linear response regime 
predicts a linear dispersion relation ujq = CsQ, with a sound velocity Cs such that 



of the order oifikp /m in the regime a < 0. The semi-classicality condition ql^cs ^ 1 then 
results in hujq ^ A. This is satisfactory, as for hujq > 2A one observes in the full RPA 
theory a coupling of the sound wave to the elementary, pair breaking excitations, which 
may distort the dispersion relation and even damp the sound wave, see e.g. [721 1731 174] 
and references therein, a phenomenon not included in superfluid hydrodynamics. 

We acknowledge useful suggestions from Felix Werner, Mattia Jona-Lasinio, Christoph 
Weiss, Alice Sinatra. 



(279) 



Aioff(r,t) = /io[p(f,t)] 



(280) 




(281) 



2 flA^Or 1 

mc, = p—[p], 



58 



YVAN CASTIN 



REFERENCES 

[1] Castin Y., "Simple theoretical tools for low dimension Bose gases", Lecture notes of the 
2003 Les Houches Spring School, Quantum Gases in Low Dimensions, Olshanii M., Perrin 
H., Pricoupenko L. Eds., J. Phys. IV France 116 (2004) 89-132. 

[2] Diu B., Guthman C, Lederer D., and Roulet B., Physique Statistique (Hermann, Paris, 
France, 1989). 

[3] Efetov K. B. and Larkin A. I., Sov. Phys. JETP 42 (1976) 390. 

[4] H.G. Vaidya and C.A. Tracy, Phys. Rev. Lett. 42, 3 (1979) and J. Math. Phys. 20, 2291 
(1979). 

[5] See the paragraph below Eq.(31) in the paper by Levitov L. S., Lee H.-W., and Lesovik 

G. B., J. Math. Phys. 37 (1996) 4845. 
[6] Wigner E. P., Phys. Rev. 40 (1932) 749. 

[7] Vignolo P., Minguzzi A., and Tosi M. P., Phys. Rev. Lett. 85 (2000) 2850. 
[8] Baranov M., private communication, October 2000. 

[9] Akdeniz Z., Vignolo P., Minguzzi A., and Tosi M. P., Phys. Rev. A 66 (2002) 055601. 
[10] Petrov D. S., Phys. Rev. Lett. 93 (2004) 143201. 

[11] Moerdijk A. J., Verhaar B. J., and Axelsson A., Phys. Rev. A 51 (1995) 4852. 

[12] Blatt J. M. and Weisskopf V. F., in Theoretical Nuclear Physics , Wiley, New York (1952). 

[13] Dalibard J., in Collisional dynamics of ultra-cold atomic gases, Proceedings of the 

International School of Physics Enrico Fermi, Course CXL: Bose - Einstein condensation 

in gases, Varenna, M. Inguscio, S. Stringari, C. Wieman edts, (1998). 
[14] Castin Y., in Coherent atomic matter waves. Lecture notes of Les Houches summer school, 

edited by Kaiser R., Westbrook C, and David F., EDP Sciences and Springer- Verlag (2001) 

1-136. 

[15] Cohen- Tannoudji C, Course at College de France, Lectures 4, 5 (1998-1999), available 

online at http:/ /www.phys. ens.fr/cours/college-de-france/1998-99/1998-99.htm 
[16] Busch T., Englert B. G., Rzazewski K., and Wilkens M., Found. Phys. 28 (1998) 549. 
[17] Idziaszek Z. and Calarco T., Phys. Rev. A 71 (2005) 050701(R). 
[18] Efimov V. N., Sov. J. Nucl. Phys. 12 (1971) 589. 
[19] Efimov V. N., Nucl. Phys. A210 (1973) 157. 

[20] Jonsell S., Heiselberg H., and Pethick C. J., Phys. Rev. Lett. 89 (2002) 250401. 
[21] Werner F. and Castin Y., Phys. Rev. Lett. 97 (2006) 150401. 
[22] Castin Y., Comptes Rendus Physique 5 (2004) 407. 
[23] Tan S., ,cond-mat/0412764, 



[24] Werner F. and Castin Y., |cond-mat/0507399| (vl) 



[25] Werner F. and Castin Y., Phys. Rev. A 74 (2006) 053604. 
[26] Petrov D. S., Phys. Rev. A 67 (2003) 010703. 

[27] Reed Michael, and Berry Simon, in Methods of modern mathematical physics 1: functional 

analysis. Academic Press (1980). 
[28] This is in a way doubly naive for a — > 0"^, as far as the pseudo-potential is concerned, 

since for a > the exact problem admits an eigenstate of energy < 3?itj/2, that — » — oo 

for a 0+. 

[29] Heiselberg H., Phys. Rev. A 63 (2001) 043606 and references therein. 

[30] Burovski E., Prokof'ev N., Svistunov B., Troyer M., Phys. Rev. Lett. 96 (2006) 160402; 

Burovski E., Prokof'ev N., Svistunov B., Troyer M., New J. Phys. 8 (2006) 153. 
[31] Bulgac A., Drut J. E., Magierski P., Phys. Rev. Lett. 96 (2006) 090404. 
[32] Lee Dean, Schaefer Thomas, Phys. Rev. C 73 (2006) 015202. 



Basic theory tools for degenerate Fermi gases 



59 



[33] Juillet O.,(cond-mat/0609063i 

[34] Mora C. and Castin Y., Phys. Rev. A 67 (2003) 053615. 

[35] Carlson J., Chang S.-Y., Pandharipande V.R. and Schmidt K. E., Phys. Rev. Lett. 91 
(2003) 050401. 

[36] Astrakharchik G. E., Boronat J., CasuUeras J., and Giorgini S., Phys. Rev. Lett. 93 (2004) 
200404. 

[37] O'Hara K. M., Hemmer S. L., Gehm M. E., Granade S. R., and Thomas J. E., Science 298 

(2002) 2179; Gehm M. E.„ Hemmer S. L., Granade S. R., O'Hara K. M., and Thomas J. 

E., Phys. Rev. A 68 (2003) 011401. 
[38] M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. Hecker Denschlag, R. Grimm, 

Phys. Rev. Lett. 92, 120401 (2004). 
[39] T. Bourdel, L. Khaykovich, J. CubizoUes, J. Zhang, F. Chevy, M. Teichmann, L. Tarruell, 

S. J. J. M. F. Kokkelmans, and C. Salomon, Phys. Rev. Lett. 93, 050401 (2004). 
[40] Greiner M., Regal C, and Jin D. S., Nature 426 (2003) 537. 

[41] Jochim S., Bartenstein M., Altmeyer A., Hendl G., Riedl S., Chin C, Denschlag J.H., and 

Grimm R., Science 302 (2003) 2101. 
[42] M.W. Zwierlein, C.A. Stan, C.H. Schunck, S.M.F. Raupach, S. Gupta, Z. Hadzibabic, and 

W. Ketterle, Phys. Rev. Lett. 91, 250401 (2003). 
[43] Partridge G. B., Li W., Kamar R. L, Liao Y. A. and Hulet R. G., Science 311 (2006) 503. 
[44] Nozieres P. and Schmitt-Rink S., J. Low Temp. Phys. 59 (1985) 195. 

[45] Randeria M., in Bose-Einstein Condensation, edited by Griffin A., Snoke D. W., and 

Stringari S. (Cambridge University Press, Cambridge) (1995) 355. 
[46] Pricoupenko L., Castin Y., Phys. Rev. A 69 (2004) 051601(R). 

[47] Cowell S., Heiselberg H., Mazets I.E., Morales J., Pandharipande V.R., Pethick C.J., Phys. 

Rev. Lett. 88 (2002) 210403. 
[48] Petrov D. S., Salomon C, and Shlyapnikov G. V., Phys. Rev. A 71 (2005) 012708; Petrov 

D. S., Salomon C, and Shlyapnikov G. V., Phys. Rev. Lett. 93 (2004) 090404. 
[49] Fulde P., Ferrell R. A., Phys. Rev. 135 (1964) A550; Larkin A.I., Ovchinnikov Y.N., Sov. 

Phys. JETP 20 (1965) 762. 
[50] Combescot R., Mora C, Eur. Phys. J. B 28 (2002) 397 and references therein. 
[51] Zwierlein M. W., Schirotzek A., Schunck C. H., and Ketterle W., Science 311 (2006) 492; 

Shin Y., Zwierlein M. W., Schunck C. H., Schirotzek A., Ketterle W., Phys. Rev. Lett. 97, 

(2006) 030401. 

[52] Partridge G. B., Li Wenhui, Liao Y. A., Hulet R. G., Haque M. and Stoof H. T. C, Phys. 

Rev. Lett. 97 (2006) 190407. 
[53] There are currently many theoretical papers on this subject. See references in e.g. [55] and 

in Chevy F., Phys. Rev. Lett. 96 (2006) 130401. 
[54] Bogoliubov N. N., Tolmachev V. V., and Shirkov D. V., New Method in the Theory of 

Superconductivity (Academy of the Sciences of the U.S.S.R., Moscow, 1958). 
[55] Anderson P. W., Phys. Rev. 112 (1958) 1900. 

[56] Jaksch D., Bruder C, Cirac J. L, Gardiner C. W., and ZoUer P., Phys. Rev. Lett. 81 
(1998) 3108. 

[57] Greiner M., Mandel O., Esslinger T., Hansch T. W., and Bloch L, Nature 415 (2002) 39. 
[58] Stoferle T., Moritz H., Giinter K., Kohl M., Esslinger T., Phys. Rev. Lett. 96 (2006) 
030401. 

[59] Chin K., Miller D. E., Liu Y., Stan C, Setiawan W., Sanner C, Xu K., Ketterle W., 

Nature 443 (2006) 961. 
[60] Bardeen J., Cooper L. N., and Schrieffer J. R., Phys. Rev. 108 (1957) 1175. 
[61] Gardiner C, ZoUer P., Quantum Noise (Springer, Heidelberg 2004). 
[62] Leggett A. J., Phys. Rev. 140 (1965) 1869. 



60 



YVAN Castin 



[63] Blaizot J. -P. and Ripka G., Quantum Theory of Finite Systems, The MIT Press 
(Cambridge, Massachusetts, 1986). 

[64] One may wonder if condition Ea. (|202|l really implies that I^'bcs) is an eigenstate of Ti. 
This is not evident a priori since the variation (5|^'bcs) does not explore the whole Hilbert 
space but only the manifold tangent to the BCS family in I'I'bcs)- One can show however 
that the answer is yes: as sketched below Ea. (|251|l . (7i — £[$0]) I'I'bcs) also belongs to this 
tangent manifold, so that one may choose (5|*I'bcs) proportional to it. 

[65] This expression coincides with Ea. (l221|l . since one may check that ([f/]~^)g- — 

[66] The reasoning was done for a fixed chemical potential and for the excitation energy of 
the grand canonical Hamiltonian H . What is the corresponding excitation energy of the 
canonical Hamiltonian for a fixed number A'^ of particles? We have a canonical energy in 
the ground state Scan = {H)^{ij.) + fiN, and in the excited state i?can = (Hyif^e) + f^eN. 
Here the chemical potentials fi and fj,e differ by 0{1/L^) for a few Bogoliubov excitations, 
so that one may expand to first order in Sfi — fie—fi. This leads to E^^^ — E^^^ = {HY{^) — 

{HY{^i) + 5^ \n + ^(-9')''(m)] +0(1/L^). Now the quantity in between brackets vanishes, 

a standard thermodynamic result. One can also use (d/dfi) |^(-?^)5^=A(fj)] ~ (^m(-^)a)a=a(m) 
and calculate the partial derivative of Ea. (|234|) with respect to fi for a fixed A. Setting A = A 
in the resulting expression leads to the standard thermodynamic result. 

[67] P. G. de Gennes, in Superconductivity of Metals and Alloys, section 4.3, Addison- Wesley 
Publishing Company (Reading, Massachusetts, fifth printing 1996). 

[68] Engelbrecht J. R., Randeria M., and Sa de Melo C, Phys. Rev. B 55 (1997) 15153. 

[69] We consider here the case of neutral particles. For charged particles the situation 
is different, because of the long range Coulomb interaction, see P. C. Martin, in 
Superconductivity, edited by R. D. Parks (Dekker, New York, 1969) Vol. 1. 

[70] Cote R. and Griffin A., Phys. Rev. B 48 (1993) 10404. 

[71] Bruun G. M. and Mottelson B. R., Phys. Rev. Lett. 87 (2001) 270403. 

[72] Minguzzi A., Ferrari G., Castin Y., Eur. Phys. J. D 17 (2001) 49. 

[73] Biichler H. P., ZoUer P., and Zwerger W., Phys. Rev. Lett. 93 (2004) 080401. 

[74] Combescot R., Kagan M. Yu., Stringari S., |cond-mat/0607493| 



[75] The RPA approach can be applied in an isotropic harmonic trap, with some non negligible 

numerical effort, see [71j . 
[76] Baranov M. A. and Petrov D. S., Phys. Rev. A 62 (2000) 041601(R). 
[77] Minguzzi A. and Tosi M. P., Phys. Rev. A 63 (2001) 023609. 
[78] Menotti C, Pedri P., Stringari S., Phys. Rev. Lett. 89 (2002) 250402. 
[79] Cozzini M., Stringari S., Phys. Rev. Lett. 91 (2003) 070401. 
[80] Tonini G., Werner F., Castin Y., Eur. Phys. J. D 39 (2006) 283. 
[81] Gor'kov L. P. and Melik-Barkhudarov T. K., Sov. Phys. JETP 13 (1961) 1018. 
[82] Zwierlein M. W., Abo-Shaeer J. R. , Schirotzek A. , Schunck C.H., Ketterle W., Nature 

435 (2005) 1047. 

[83] Messiah A., in Quantum Mechanics, vol.11 (North Holland, 1961); Migdal A.B., in 
Qualitative methods in Quantum Theory, (W.A. Benjamin, Massachusetts, 1977) 155. 

[84] D. R. Meacher, D. Boiron, H. Metcalf, C. Salomon, and G. Grynberg, Phys. Rev. A 50 
(1994) R1992 - R1994 . 

[85] Stamper-Kurn D. M., Chikkatur A. P., Gorlitz A., Inouye S., Gupta S., Pritchard D. E. 
and Ketterle W., Phys. Rev. Lett. 83 (1999) 2876. 



